Adjoint transport for dose in beam angle optimization for external beam radiation therapy

ABSTRACT

A method of beam angle optimization for an IMRT radiotherapy treatment includes providing a patient model having one or more regions of interest (ROIs), defining a delivery coordinate space (DCS), for each ROI, solving an adjoint transport to obtain an adjoint solution field from the ROI, for each vertex in the DCS, evaluating an adjoint photon fluence by performing ray tracing of the adjoint solution field, evaluating a dose of the ROI using the adjoint photon fluence, for each vertex in the DCS, evaluating a respective beam&#39;s eye view (BEV) score of each pixel of a BEV plane using the doses of the one or more ROIs, determining one or more BEV regions in the BEV plane based on the BEV scores, determining a BEV region connectivity manifold based on the BEV regions, and determining a set of IMRT fields based on the BEV region connectivity manifold.

CROSS-REFERENCE TO RELATED APPLICATIONS

The present application is a non-provisional application of and claimsthe benefit and priority under 35 U.S.C. 119(e) of U.S. ProvisionalApplication No. 62/738,741, filed Sep. 28, 2018, entitled “ADJOINTTRANSPORT FOR DOSE IN TREATMENT TRAJECTORY OPTIMIZATION AND BEAM ANGLEOPTIMIZATION FOR EXTERNAL BEAM RADIATION THERAPY,” the entire content ofwhich is incorporated herein by reference for all purposes.

The following three U.S. patent applications (including this one) arebeing filed concurrently, and the entire disclosures of the otherapplications are incorporated by reference into this application for allpurposes:

U.S. application Ser. No. 16/586,654, filed Sep. 27, 2019, entitled“ADJOINT TRANSPORT FOR DOSE IN TREATMENT TRAJECTORY OPTIMIZATION FOREXTERNAL BEAM RADIATION THERAPY”;

U.S. application Ser. No. 16/586,661, filed Sep. 27, 2019, entitled“ADJOINT TRANSPORT FOR DOSE IN BEAM ANGLE OPTIMIZATION FOR EXTERNAL BEAMRADIATION THERAPY”; and

U.S. application Ser. No. 16/586,666, filed Sep. 27, 2019, entitled“HYBRID TRAJECTORY AND BEAM ANGLE OPTIMIZATION FOR EXTERNAL BEAMRADIATION THERAPY”.

BACKGROUND

Significant developments have been made in inverse treatment planning ofexternal beam radiation therapies using, e.g., IMRT and VMAT treatmentmodalities. As both plan quality requirements and requirements on theclinics' patient throughput increase, the role of automation and ahigher degree of personalization of treatment plans become increasinglyimportant. Prior to inverse-optimization of the MLC leaf sequence anddose rates in a radiation treatment plan, treatment fields, such as VMATtrajectories and IMRT fields, may need to be determined. Each treatmentfield may be associated with a beam energy (e.g., 6 MeV or 15 MeV) andprofile (e.g., flat or flattening-filter free). The choice of beamenergy and treatment geometry (e.g., isocenter(s), starting and stoppinggantry angles, and collimator angle(s) of VMAT trajectories, or gantryangle, couch angle, collimator angle, and jaw positions of IMRT fields)may be dictated by a clinical protocol for a given treatment site. Sucha protocol may be sub-optimal considering the variability in patientanatomy and in clinical goals. Therefore, there is a need for improvedmethods of optimizing treatment geometries.

SUMMARY

According to some embodiments, a method of trajectory optimization forradiotherapy treatment includes providing a patient model that includesone or more regions of interest (ROIs) for the radiotherapy treatment,and defining a delivery coordinate space (DC S) having a set ofcandidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane. The method further includes, foreach respective ROI of the one or more ROIs, solving an adjointtransport to obtain an adjoint solution field from the respective ROI;and for each respective candidate vertex in the DCS, and for eachrespective pixel of the respective BEV plane defined by the respectivecandidate vertex, evaluating an adjoint photon fluence originating froma respective beamlet incident from the respective candidate vertex andpassing through the respective pixel by performing ray tracing of theadjoint solution field; and evaluating a respective dose of therespective ROI from the respective beamlet using the adjoint photonfluence. The method further includes for each respective candidatevertex in the DCS, and for each respective pixel of the respective BEVplane defined by the respective candidate vertex, evaluating arespective BEV score of the respective pixel using the doses of the oneor more ROIs evaluated for the respective beamlet incident from therespective candidate vertex and passing through the respective pixel;and determining one or more BEV regions in the respective BEV planebased on the BEV scores of the pixels of the BEV planes corresponding tothe set of candidate vertices in the DCS. The method further includesdetermining a BEV region connectivity manifold based on the BEV regionsof the BEV planes of the set of candidate vertices in the DCS. The BEVregion connectivity manifold represents connections between contiguousBEV regions between adjacent vertices. The method further includesdetermining one or more optimal treatment trajectories based on the BEVregion connectivity manifold.

According to some embodiments, a method of trajectory optimization forradiotherapy treatment includes providing a patient model that has oneor more regions of interest (ROIs) for the radiotherapy treatment, anddefining a delivery coordinate space (DCS) having a set of candidatevertices. Each respective candidate vertex defines a respective beam'seye view (BEV) plane. The method further includes identifying aplurality of candidate energy modes for the radiotherapy treatment; andfor each respective energy mode of the plurality of candidate energymodes, for each respective BEV plane of a respective candidate vertex inthe DCS, and for each respective ROI of the one or more ROIs, evaluatinga respective dose of the respective ROI from a respective beamletincident from the respective candidate vertex and passing through arespective pixel of the respective BEV plane using transport solutionsfor the respective energy mode. The method further includes evaluating arespective BEV score of the respective pixel using the doses of the oneor more ROIs evaluated for the respective beamlet; determining one ormore BEV regions in the respective BEV plane for the respective energymode based on the BEV scores of the pixels of the BEV planescorresponding to the set of candidate vertices in the DCS; anddetermining a respective BEV region connectivity manifold for therespective energy mode based on the BEV regions of the BEV planes of theset of candidate vertices in the DCS. The respective BEV regionconnectivity manifold represents connections between contiguous BEVregions between adjacent vertices. The method further includesdetermining a plurality of candidate sets of optimal treatmenttrajectories by: for each respective energy mode of the plurality ofcandidate energy modes, determining a respective candidate set ofoptimal treatment trajectories based on the respective BEV regionconnectivity manifold for the respective energy mode. The method furtherincludes selecting one of the plurality of candidate sets of optimaltreatment trajectories as a final set of optimal treatment trajectoriesbased on an objective function. The final set of optimal treatmenttrajectories corresponds to an optimal energy mode among the pluralityof candidate energy modes.

According to some embodiments, a method of trajectory optimization forradiotherapy treatment includes providing a patient model including oneor more regions of interest (ROIs) for the radiotherapy treatment, anddefining a delivery coordinate space (DCS) having a set of candidatevertices. Each respective candidate vertex defines a respective beam'seye view (BEV) plane. The method further includes identifying a firstenergy mode and a second energy mode for the radiotherapy treatment; forthe first energy mode, determining a first BEV region connectivitymanifold by evaluating doses of the one or more ROIs for each respectiveBEV plane corresponding to each respective candidate vertex of the DCSusing transport solutions for the first energy mode; and for the secondenergy mode, determining a second BEV region connectivity manifold byevaluating doses of the one or more ROIs for each respective BEV planecorresponding to each respective candidate vertex of the DCS usingtransport solutions for the second energy mode. The method furtherincludes determining a set of optimal treatment trajectories based onthe first BEV region connectivity manifold and the second BEV regionconnectivity manifold.

According to some embodiments, a method of beam angle optimization foran IMRT radiotherapy treatment include providing a patient model thathas one or more regions of interest (ROIs) for the IMRT radiotherapytreatment, and defining a delivery coordinate space (DCS) having a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane. The method further includes, foreach respective ROI of the one or more ROIs, solving an adjointtransport to obtain an adjoint solution field from the respective ROI;and for each respective candidate vertex in the DCS, for each respectivepixel of the respective BEV plane defined by the respective candidatevertex, evaluating an adjoint photon fluence originating from arespective beamlet incident from the respective candidate vertex andpassing through the respective pixel by performing ray tracing of theadjoint solution field; evaluating a respective dose of the respectiveROI from the respective beamlet using the adjoint photon fluence. Themethod further includes, for each respective candidate vertex in theDCS, and for each respective pixel of the respective BEV plane definedby the respective candidate vertex, evaluating a respective BEV score ofthe respective pixel using the doses of the one or more ROIs evaluatedfor the respective beamlet incident from the respective candidate vertexand passing through the respective pixel; and determining one or moreBEV regions in the respective BEV plane based on the BEV scores of thepixels of the BEV planes corresponding to the set of candidate verticesin the DCS. The method further includes determining a BEV regionconnectivity manifold based on the BEV regions of the BEV planes of theset of candidate vertices in the DCS. The BEV region connectivitymanifold represents connections between contiguous BEV regions betweenadjacent vertices. The method further includes determining a set of IMRTfields based on the BEV region connectivity manifold. Each respectiveIMRT field of the set of IMRT fields defines a beam angle correspondingto a respective vertex in the DCS.

According to some embodiments, a method of beam angle optimization foran IMRT radiotherapy treatment includes providing a patient model havingone or more regions of interest (ROIs) for the IMRT radiotherapytreatment, and defining a delivery coordinate space (DCS) having a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane. The method further includesidentifying a plurality of candidate energy modes for the IMRTradiotherapy treatment; and for each respective energy mode of theplurality of candidate energy modes, for each respective BEV plane of arespective candidate vertex in the DCS, and for each respective ROI ofthe one or more ROIs: evaluating a respective dose of the respective ROIfrom a respective beamlet incident from the respective candidate vertexand passing through a respective pixel of the respective BEV plane usingtransport solutions for the respective energy mode; evaluating arespective BEV score of the respective pixel using the doses of the oneor more ROIs evaluated for the respective beamlet; and determining oneor more BEV regions in the respective BEV plane for the respectiveenergy mode based on the BEV scores of the pixels of the BEV planescorresponding to the set of candidate vertices in the DCS. The methodfurther includes determining a respective BEV region connectivitymanifold for the respective energy mode based on the BEV regions of theBEV planes of the set of candidate vertices in the DCS. The respectiveBEV region connectivity manifold represents connections betweencontiguous BEV regions between adjacent vertices. The method furtherincludes determining a plurality of candidate sets of IMRT fields by,for each respective energy mode of the plurality of candidate energymodes, determining a respective candidate set of IMRT fields based onthe respective BEV region connectivity manifold for the respectiveenergy mode, each respective IMRT field defining a beam anglecorresponding to a respective vertex in the DCS. The method furtherincludes selecting one of the plurality of candidate sets of IMRT fieldsas an optimal set of IMRT fields based on an objective function. Theoptimal set of IMRT fields corresponds to an optimal energy mode amongthe plurality of candidate energy modes.

According to some embodiments, a method of beam angle optimization in anIMRT radiotherapy treatment includes providing a patient model includingone or more regions of interest (ROIs) for the IMRT radiotherapytreatment, and defining a delivery coordinate space (DCS) having a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane. The method further includesidentifying a first energy mode and a second energy mode for the IMRTradiotherapy treatment; for the first energy mode, determining a firstBEV region connectivity manifold by evaluating doses of the one or moreROIs for each respective BEV plane corresponding to each respectivecandidate vertex of the DCS using transport solutions for the firstenergy mode; and for the second energy mode, determining a second BEVregion connectivity manifold by evaluating doses of the one or more ROIsfor each respective BEV plane corresponding to each respective candidatevertex of the DCS using transport solutions for the second energy mode.The method further includes determining a set of IMRT fields based onthe first BEV region connectivity manifold for the first energy mode andthe second BEV region connectivity manifold for the second energy mode.Each respective IMRT field of the set of IMRT fields defines a beamangle corresponding to a respective vertex in the DCS.

According to some embodiments, a method of determining treatmentgeometries for a radiotherapy treatment includes providing a patientmodel having one or more regions of interest (ROIs) for the radiotherapytreatment, and defining a delivery coordinate space (DCS) having a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane. The method further includes, foreach respective BEV plane of a respective candidate vertex in the DCS,and for each respective ROI of the one or more ROIs, evaluating arespective dose of the respective ROI from a respective beamlet incidentfrom the respective candidate vertex and passing through a respectivepixel of the respective BEV plane using transport solutions of therespective beamlet; evaluating a respective BEV score of the respectivepixel using the doses of the one or more ROIs evaluated for therespective beamlet; and determining one or more BEV regions in therespective BEV plane based on the BEV scores of the pixels of the BEVplanes corresponding to the set of candidate vertices in the DCS. Themethod further includes determining a BEV region connectivity manifoldbased on the BEV regions of the BEV planes of the set of candidatevertices in the DCS. The BEV region connectivity manifold representsconnections between contiguous BEV regions between adjacent vertices.The method further includes determining a set of treatment trajectoriesbased on the BEV region connectivity manifold. Each treatment trajectorydefines a respective path through a respective set of vertices in theDCS. The method further includes determining one or more IMRT fields,each respective IMRT field defines a respective direction of incidencecorresponding to a respective vertex in the DCS.

According to some embodiments, a method of determining treatmentgeometries for a radiotherapy treatment includes providing a patientmodel having one or more regions of interest (ROIs) for the radiotherapytreatment, and defining a first delivery coordinate space and a seconddelivery coordinate space. The first delivery coordinate space has afirst set of candidate vertices. The second delivery coordinate spacehas a second set of candidate vertices. Each vertex of the first set ofcandidate vertices or the second set of candidate vertices defines arespective beam's eye view (BEV) plane. The method further includesdetermining a first beam's eye view (BEV) region connectivity manifoldby evaluating doses of the one or more ROIs for each respective BEVplane corresponding to each respective vertex of the first set ofcandidate vertices of the first delivery coordinate space, anddetermining a first set of treatment trajectories based on the first BEVregion connectivity manifold. Each treatment trajectory of the first setof treatment trajectories defines a respective path through a respectiveset of vertices in the first delivery coordinate space. The methodfurther includes determining a first set of IMRT fields. Each of thefirst set of IMRT fields corresponds to a respective vertex in thesecond delivery coordinate space.

According to some embodiments, a method of determining treatmentgeometries for a radiotherapy treatment includes providing a patientmodel having one or more regions of interest (ROIs) for the radiotherapytreatment, and defining a delivery coordinate space (DCS) having a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane. The method further includesidentifying a first energy mode and a second energy mode for theradiotherapy treatment; for the first energy mode, determining a firstBEV region connectivity manifold by evaluating doses of the one or moreROIs for each respective BEV plane corresponding to each respectivecandidate vertex of the DCS using transport solutions for the firstenergy mode; for the second energy mode, determining a second BEV regionconnectivity manifold by evaluating doses of the one or more ROIs foreach respective BEV plane corresponding to each respective candidatevertex of the DCS using transport solutions for the second energy mode;determining a set of optimal treatment trajectories based on the firstBEV region connectivity manifold and the second BEV region connectivitymanifold; and determining a set of IMRT fields. Each respective IMRTfield of the set of IMRT fields corresponds to a respective vertex inthe delivery coordinate space.

These and other embodiments of the disclosure are described in detailbelow. For example, other embodiments are directed to systems, devices,and computer readable media associated with methods described herein.

A better understanding of the nature and advantages of embodiments ofthe present disclosure may be gained with reference to the followingdetailed description and the accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed incolor. Copies of this patent or patent application publication withcolor drawing(s) will be provided by the Office upon request and paymentof the necessary fee.

FIG. 1 is a schematic perspective view of a radiation treatment system.

FIG. 2 is a schematic side view of a radiation treatment system.

FIG. 3 shows schematically a photon collimation system in a radiationtreatment system.

FIG. 4 shows an exemplary multi-leaf collimator (MLC) plane.

FIG. 5 shows a block diagram of an external-beam radiation treatmentsystem of FIGS. 1 and 2 .

FIG. 6 shows a representation of a delivery coordinate space (DCS)truncated to avoid collisions, where each point (θ_(gantry), θ_(couch))is transformed to the physical three-dimensional (3D) location of atreatment head.

FIG. 7 shows examples of a differential beam's eye view (BEV) fiberbundle section histogram and an integrated BEV fiber bundle sectionhistogram according to some embodiments.

FIGS. 8A and 8B show exemplary images of BEV scores for two adjacentdelivery coordinate vertices for a chest wall treatment, respectfully,according to some embodiments.

FIG. 9 illustrates an exemplary quadrilateralized spherical cube thatmay be used for evaluating the angular flux according to someembodiments.

FIG. 10A illustrates boundary nodes in a one-dimensional (1D) deliverycoordinate space according to some embodiments.

FIG. 10B illustrates boundary nodes in a two-dimensional (2D) deliverycoordinate space according to some embodiments.

FIG. 11 shows a flowchart illustrating a method of trajectoryoptimization according to some embodiments.

FIG. 12 shows a flowchart illustrating a method of beam angleoptimization according to some embodiments.

FIGS. 13A and 13B show exemplary BEV regions for a given direction ofincidence for two different energy modes, respectively, according tosome embodiments.

FIG. 14 is a flowchart illustrating a method of trajectory optimizationconsidering multiple candidate energy modes according to someembodiments.

FIG. 15 is a flowchart illustrating a method of trajectory optimizationconsidering multiple candidate energy modes according to someembodiments.

FIG. 16 is a flowchart illustrating a method of beam angle optimizationconsidering multiple candidate energy modes according to someembodiments.

FIG. 17 is a flowchart illustrating a method of beam angle optimizationconsidering multiple candidate energy modes according to someembodiments.

FIG. 18A illustrates a configuration of a C-arm linear acceleratorradiotherapy treatment system according to some embodiments.

FIG. 18B shows a view of a patient model from the top of the patient'shead and two exemplary treatment trajectories according to someembodiments.

FIG. 18C shows a view of a patient model from the top of the patient'shead and two exemplary treatment trajectories according to someembodiments.

FIG. 19 shows an example of a more general treatment trajectorysuperimposed on a patient model according to some embodiments.

FIG. 20 is a flowchart illustrating a method of hybrid VMAT-IMRTtreatment geometry optimization according to some embodiments.

FIG. 21A illustrates a first delivery coordinate space in a firstconfiguration of a C-arm linear accelerator radiotherapy treatmentsystem according to some embodiments.

FIG. 21B shows a treatment trajectory in the first delivery coordinateillustrated in FIG. 21A according to some embodiments.

FIG. 21C illustrates a second delivery coordinate space in a secondconfiguration of a C-arm linear accelerator radiotherapy treatmentsystem according to some embodiments.

FIG. 21D shows a treatment trajectory in the second delivery coordinateillustrated in FIG. 21C according to some embodiments.

FIG. 22 shows two treatment trajectories in a third delivery coordinatespace and a fourth delivery coordinate space, respectively.

FIG. 23 is a flowchart illustrating a method of hybrid VMAT-IMRTtreatment geometry optimization according to some embodiments.

FIG. 24 is a flowchart illustrating a method of hybrid VMAT-IMRTtreatment geometry optimization according to some embodiments.

FIG. 25 shows a block diagram of an example computer system usable withsystem and methods according to embodiments.

DETAILED DESCRIPTION

The present disclosure relates generally to treatment planning forradiotherapy treatment using external-beam radiotherapy treatmentsystems, and is more particularly directed to optimizing trajectoriesand field geometries in a radiation treatment plan. Beam's eye view(BEV) regions and BEV region connectivity manifold may be determined byevaluating dose response of each region of interest for each vertex in adelivery coordinate space (DCS). The information contained in the BEVregions and BEV region connectivity manifold may be used to generateoptimized trajectories or optimized field geometries in a radiationtreatment plan.

The detailed description of the various embodiments is organized asfollows. Section I discusses some exemplary radiation treatment systems.Section II discusses radiation treatment planning. Section III discussesbeam's eye view (BEV) sectioning for an approach to treatment trajectoryoptimization in radiotherapy (referred to herein as TORUS). The TORUSapproach may require computationally costly and time-consuming doseevaluations. Section IV discusses two approaches to dose evaluation: theforward transport approach and the adjoint transport approach. Theadjoint transport approach may have the advantage of morecomputationally efficient dose evaluation, especially for optimizationsthat consider multiple candidate energy modes. Section V discussestrajectory optimization using adjoint transport for dose. Section VIdiscusses beam angle optimization in IMRT radiotherapy treatment usingadjoint transport for dose. Section VII discusses trajectoryoptimization considering multiple candidate energy modes. Section VIIIdiscusses beam angle optimization considering multiple candidate energymodes. Section IX discusses applications of TORUS and BEV sectioningbased methodologies to hybrid trajectory and beam angle optimization.

I. Treatment System

FIGS. 1 and 2 depict a radiation treatment system of the type that maybe used in connection with the present invention. While the inventionsherein will be described with reference to C-arm gantries, it will beappreciated by those skilled in the art that the same principles applyto other rotatable gantries such as ring gantries. Referring to FIG. 1 ,a perspective view of radiation treatment system (in this case a linearaccelerator) is shown. Typically, such a system is capable of generatingeither an electron (particle) beam or an x-ray (photon) beam for use inthe radiotherapy treatment of patients on a treatment couch 35. Otherradiation treatment systems are capable of generating heavy ionparticles such as protons. For purposes of the present discussion, onlyx-ray irradiation will be discussed. However, it will be appreciated bythose skilled in the art that the same principles apply to othersystems.

Stand 10 supports a rotatable gantry 20 with a treatment head 30. Nextto stand 10 there is arranged a control unit (not shown) that includescontrol circuitry for controlling the different modes of operation ofthe accelerator. A high voltage source is provided within the stand orin the gantry, to supply voltage to an electron gun (not shown)positioned on an accelerator guide located in the gantry 20. Electronsare emitted from the electron gun into the guide (not shown) where theyare accelerated. A source supplies RF (microwave) power for thegeneration of an electric field within the waveguide. The electronsemitted from the electron gun are accelerated in the waveguide by theelectric field, and exit the waveguide as a high energy electron beam,typically at megavoltage energies. The electron beam then strikes asuitable metal target, emitting high energy x-rays in the forwarddirection.

Referring now to FIG. 2 , a somewhat more detailed side view of aradiation treatment system of the type that may be used in connectionwith the present invention is shown. A patient P is shown lying on thetreatment couch 35. X-rays formed as described above are emitted fromthe target in the treatment head 30 in a divergent beam 104. Typically,a patient plane 116, which is perpendicular to the page in FIG. 2 , ispositioned about one meter from the x-ray source or target, and the axisof the gantry 20 is located on the plane 116, such that the distancebetween the target and the isocenter 178 remains constant when thegantry 20 is rotated. The isocenter 178 is at the intersection betweenthe patient plane 116 and the central axis of beam 122. A treatmentvolume to be irradiated is located about the isocenter 178.

FIG. 3 shows schematically a photon collimation system 300 with upperjaws 310 (i.e., the Y1 and Y2 jaws; the Y1 jaw is omitted for clarity),lower jaws 320 (i.e., the X1 and X2 jaws), and a multi-leaf collimator(MLC) 330. The field dimensions in the plane 340 at the isocenter 178are indicated. The upper jaws 310, the lower jaws 320, and the leaves332 of the MLC 330 comprise an x-ray blocking material, and arepositioned in the head 30 to define the width of the x-ray beam at thepatient plane. Typically, the jaws 310 and 320 are moveable and, whenfully open, define a maximum beam of about 40 cm×40 cm at the patientplane 116. The MLC 330 is positioned at the exit of the head 30, tofurther shape the x-ray beam. Since its introduction in 1990 the MLC hasbecome a standard feature of most radiation treatment systems. Anexample of a current MLC sold by the assignee of the present inventionuse up to 120 individually controllable leaves, typically thin slices oftungsten, that can be moved into or out of the x-ray beam under thecontrol of system software.

FIG. 4 shows an exemplary MLC plane having a plurality of leaves 332,arranged in opposing pairs, and an aperture 415 created by selected leafmovements. Radiation passes through and is shaped by the aperture 415.Thus, the MLC can be used to collimate the x-rays to provide conformaltreatment of tumors from various angles (“3D conformal”) as well asintensity modulated radiotherapy (“IMRT”), whereby different radiationdoses are delivered to different portions of the treatment area. Thetreatment volume, i.e., the irradiated volume proximate to the isocenter178 in the path of the x-ray beam, is defined by the jaws 310 and 320,the leaf sequences of the MLC 330, and the collimator angle, i.e., theangle at which the MLC 330 sits in the head 30. Some external radiationtreatment systems may include multiple layers of MLCs. The multiplelayers of MLCs may be positioned at different planes and at differentcollimator angles.

FIG. 5 shows a block diagram of an external-beam radiation treatmentsystem 500 of FIGS. 1 and 2 . The radiation treatment system 500includes a beam source 510, a beam aperture 520, a gantry 530, and acouch 540. The beam source 510 is configured to generate a beam oftherapeutic radiation. This beam of radiation may include x-rays,particles, and the like. The beam aperture 520 includes an adjustablemulti-leave collimator (MLC) 522 for spatially filtering the radiationbeam. The couch 540 is configured to support and position a patient. Thecouch 540 may have six degrees of freedom, namely the translationaloffsets X, Y, and Z, and the rotation, pitch, and yaw.

The gantry 530 that circles about the couch 540 houses the beam source510 and the beam aperture 520. The beam source 510 is optionallyconfigured to generate imaging radiation as well as therapeuticradiation. The radiation treatment system 500 may further include animage acquisition system 550 that comprises one or more imagingdetectors mounted to the gantry 530.

The radiation treatment system 500 further includes a control circuitry560 for controlling the operation of the beam source 510, the beamaperture 520, the gantry 530, the couch 540, and the image acquisitionsystem 550. The control circuitry 560 may include hardware, software,and memory for controlling the operation of these various components ofthe radiation treatment system 500. The control circuitry 560 cancomprise a fixed-purpose hard-wired platform or can comprise a partiallyor wholly-programmable platform. The control circuitry 560 is configuredto carry out one or more steps, actions, and other functions describedherein. In some embodiments, the control circuitry 560 may include amemory for receiving and storing a radiation treatment plan that definesthe control points of one or more treatment fields. The controlcircuitry 560 may then send control signals to the various components ofthe radiation treatment system 500, such as the beam source 510, thebeam aperture 520, the gantry 530, and the couch 540, to execute theradiation treatment plan. In some embodiments, the control circuitry 560may include an optimization engine 562 configured for determining aradiation treatment plan. In some other embodiments, the controlcircuitry 560 may not include an optimization engine. In those cases, aradiation treatment plan may be determined by an optimization engine ina separate computer system, and the radiation treatment plan is thentransmitted to the control circuitry 560 of the radiation treatmentsystem 500 for execution.

II. Radiation Treatment Planning

Radiation therapy is generally implemented in accordance with aradiation treatment plan that typically takes into account the desireddose of radiation that is prescribed to be delivered to the tumor, aswell as other constraints on dose in the surrounding tissues that dependon the tissue type. Various techniques for developing radiationtreatment plans may be used. Preferably, the computer system used todevelop the radiation treatment plan provides an output that can be usedto control the radiation treatment system, including the control pointsand the MLC leaf movements. Typically, the desired dose prescribed in aradiation treatment plan is delivered over several sessions, calledfractions.

Several techniques have been developed to create radiation treatmentplans for IMRT or conformal radiation therapy. Generally, thesetechniques are directed to solving the “inverse” problem of determiningthe optimal combination of angles, radiation doses and MLC leafmovements to deliver the desired total radiation dose to the targetwhile minimizing irradiation of healthy tissue. This inverse problem iseven more complex for developing arc therapy plans, such as volumetricmodulated arc therapy (VMAT), where the one or more external treatmentcoordinates, such as the isocenter location, gantry angle, couch angles,and couch offsets, are in motion while irradiating the target volume.Heretofore, radiation oncologists or other medical professionals, suchas medical physicists and dosimetrists, have used one of the availablealgorithms to develop and optimize a radiation treatment plan.

Typically, such planning starts with volumetric information about thetarget tumor and about any nearby tissue structures. For example, suchinformation may comprise a map of the planning target volume (“PTV”),such as a prostate tumor, which is prescribed by the physician toreceive a certain therapeutic radiation dose with allowable tolerances.Volumetric information about nearby tissues may include for example,maps of the patient's bladder, spinal cord and rectum, each of which maybe deemed an organ at risk (OAR) that can only receive a much lower,maximum prescribed amount of radiation. This volumetric informationalong with the prescribed dose limits and similar objectives set by themedical professionals are the basis for calculating an optimized dosedistribution, also referred to as fluence map, which in turn is thebasis for determining a radiation treatment plan. The volumetricinformation may, for example, be reduced to an objective function or asingle figure of merit that accounts for the relative importance ofvarious tradeoffs inherent in a radiation treatment plan, along withconstraints that must be met for the radiation treatment plan to bemedically acceptable or physically possible.

III. Beam's Eye View (BEV) Sectioning in Treatment Geometry Optimization

State-of-the-art techniques for optimizing treatment trajectories (e.g.,VMAT trajectories) in external beam radiotherapy involve dosimetriccharacterization of candidate directions of incidence. (See e.g.,Christopher Barry Locke and Karl Kenneth Bush, Trajectory OptimizationIn Radiotherapy Using Sectioning (referred to herein as TORUS), MedicalPhysics, 2017, and U.S. patent application Ser. Nos. 16/235,205 and16/235,211) The goal of the optimizations may be to ascertain whichdirections of incidence in the permissible delivery coordinate space aremore suited for treating the patient, considering the dose response ofboth planning target volumes (PTVs) and organs at risk (OARs) within thepatient.

A. Delivery Coordinate Space (DCS)

The delivery coordinate space (DCS) is a set of all allowablecoordinates that parameterize the delivery device's configuration,truncated to avoid collisions (e.g., machine-to-machine collisions andmachine-to-patient collisions). For a C-arm linear accelerator withfixed isocenter, points in a delivery coordinate space may be defined astuples of the form (θ_(gantry), θ_(couch)), where θ_(gantry) is thegantry angle, and θ_(couch) is the couch angle. The DCS may bediscretized into a 2D mesh defined by a set of vertices (e.g., eachvertex having associated gantry angle and couch angle values), edges,and triangle faces. Thus, the DCS

may be represented by a simplicial complex (mesh) defined by:

vertices: tuples of (θ_(gantry), θ_(couch)) angles;

edges: ordered pair of start and end vertices; and

triangles: ordered list of 3 vertex indices.

FIG. 6 shows a representation of a DCS as a 3D mesh, where each point610 (θ_(gantry), θ_(couch)) is transformed to the physicalthree-dimensional (3D) location of the treatment head, and the space istruncated to avoid collisions.

B. BEV Sectioning

Christopher Barry Locke and Karl Kenneth Bush, Trajectory OptimizationIn Radiotherapy Using Sectioning (referred to herein as TORUS), MedicalPhysics, 2017 discusses a method of trajectory optimization inradiotherapy using sectioning. Dosimetry experience in radiationtreatment planning has shown that BEV offers a valuable tool indetermining the geometrical setup for both dynamic gantry treatment(e.g., VMAT) and static gantry treatment (e.g., IMRT).

For a source position r_(v) ^(s) corresponding to vertex v in the 3DDCS, a BEV plane (also referred to as an isocenter plane) may be definedas a plane perpendicular to the vector r_(v) ^(s)−r_(v) ^(ISO) andincluding the isocenter r_(v) ^(ISO). A BEV plane may be discretizedinto a 2D array of N_(x)×N_(y) pixels, with each pixel on this 2D gridrepresenting a single beamlet.

To probe all possible beamlets, the intensity of each beamlet may be setto unity, and a 3D dose response to each beamlet may be evaluated. The3D dose may be processed to determine dose statistics to each region ofinterest (ROI) for each beamlet. The ROIs may include, for example,planning target volumes (PTVs) and organs at risk (OARs). If the dosefor a beamlet at pixel (n_(x), n_(y)) in the BEV plane given by deliverycoordinate vertex n_(v) to the ROI with index n_(ROI) at position (x, y,z) is given by D_(n) _(v) _(,n) _(x) _(,n) _(y) _(,n) _(ROI) (x, y, z),then the volume integrated dose for the ROI with index n_(ROI) may beevaluated as,

_(n) _(v) _(,n) _(x) _(,n) _(y) _(,n) _(ROI) =∫dxdydzD _(n) _(v) _(,n)_(x) _(,n) _(y) _(,n) _(ROI) (x,y,z).  (1)

If there are N_(ROI) regions of interest, then the BEV dose bundlesection

is a 4D array of size (

, N_(x), N_(y), N_(ROI)) containing the volume integrated dose valuesfor each ROI from each beamlet.

In some embodiments, a BEV score bundle section

may be defined as a contraction of the 4D BEV dose bundle section

into a 3D matrix of size (

, N_(x), N_(y)), where the values are a measure of the “goodness” ofbeamlets based on the ROI dosimetrics. In some embodiments, the“goodness” score may be evaluated as a linear combination of the dosesto each ROI,

$\begin{matrix}{\mathcal{S}_{n_{v},n_{x},n_{y}} = {{\sum}_{n_{ROI} = 1}^{N_{ROI}}w_{n_{ROI}}{\mathcal{D}_{n_{v},n_{x},n_{y},n_{ROI}}.}}} & (2)\end{matrix}$

The coefficients w_(n) _(ROI) may be set by users. For example, thevalues of the coefficients may be set to −0.2 for body, between −1 and−10 for an OAR (e.g., critical organs may be given more negativeweights), and between zero and +1 for a PTV.

C. BEV Regions and BEV Region Connectivity Manifold

A BEV region connectivity manifold may be constructed in two steps.First, information contained in the BEV score bundle sections may beconsidered and a binary selection procedure is applied to determine if agiven pixel (beamlet) is a “good” or “bad” candidate for treatment. Foreach BEV plane, a set of “good” beamlets form regions. Each regionincludes a set of contiguous pixels in the BEV plane, and representpotential open aperture candidates for use in the optimization. Next,how the regions connect to other regions in neighboring vertices may bedetermined. The resulting structure, comprised of regions and theirconnections, forms a BEV region connectivity manifold.

1. BEV Region Score

In some embodiments, a beamlet may be deemed a “good” candidate if itintersects a PTV and its score S is above a certain threshold

_(threshold). Choosing an appropriate threshold may be a non-trivialtask and can be case specific. For example, beamlets treating asuperficial target with very little body or OAR in the way (e.g. a pronebreast irradiation) may have a different threshold than a deep seatedtarget (e.g., in a prostate treatment) in which the best possible planmay still treat through healthy tissue to a substantial depth.

According to some embodiments, a region score

_(n) _(v) _(,n) _(x) _(,n) _(y) ∈(−∞, 1] may be used to define regions,where the potentially “good” beamlets are defined to have

>0,

$\begin{matrix}{\mathcal{R}_{n_{v},n_{x},n_{y}} = \left\{ {\begin{matrix}{- \infty} & {{if}{\nexists{n_{ROI} \in N_{PTV} \ni {\mathcal{D}_{n_{v},n_{x},n_{y},n_{ROI}} > 0}}}} \\\frac{\mathcal{S}_{n_{v},n_{x},n_{y}} - \mathcal{S}_{threshold}}{{\max(\mathcal{S})} - \mathcal{S}_{threshold}} & {{othe}rwise}\end{matrix},} \right.} & (3)\end{matrix}$

where N_(PTV) is the set of ROI indices for PTV regions of interest.Thus, a beamlet may be deemed a “good” candidate if it intersects thePTV and its score is above some threshold

_(threshold). The region score

may classify beamlets into regions, and may also act as a normalizedscore for the goodness of beamlets (e.g., the maximum region score beingunity).

2. Score Threshold Determination

According to some embodiments, the score threshold

_(threshold) may be automatically determined using histograms of the BEVfiber bundle sections, in the spirit of dose-volume histograms (DVHs).Given a section

, where n is an index in some set

, and a subset of indices under consideration

, the associated BEV fiber bundle section histogram may be defined asfollows:

-   -   Determine the maximum and minimum values of the set {        |n∈        }, F_(max) and F_(min).    -   Create N_(bins) that range from F_(min) to F_(max), and        initialize each to 0. These may be referred to as differential        histogram bins and are denoted        for each bin index n_(bin).    -   Loop through each n∈        and increment the bin that corresponds to the value        by 1.    -   The integrated histogram bins may be defined as

$B_{n_{bin}}^{\mathcal{F}} = {{\sum}_{n = n_{bin}}^{N_{bins} - 1}{{\partial B_{n}^{\mathcal{F}}}.}}$

-   -   Normalize the differential and integrated histogram bins to have        maximum value of 1.

FIG. 7 shows examples of a differential BEV fiber bundle sectionhistogram 710 (dashed line) and an integrated BEV fiber bundle sectionhistogram 720 (solid line) according to some embodiments. The dottedlines 730 meet at the integrated histogram curve 720 and may have avertical height to horizontal length ratio of 0.7. Their intersectionwith the abscissa may give the score threshold

_(threshold).

Using this definition, one may define the BEV PTV dose histogram fromthe BEV fiber bundle section

_(n) _(v) _(,n) _(x) _(,n) _(y) ^(PTV)=Σ_(n) _(ROI) _(∈N) _(PTV)

_(n) _(v) _(,n) _(x) _(,n) _(y) _(,n) _(ROI) (the BEV fiber bundlesection representing the total PTV dose only). These histograms may bedenoted as

and

. In some embodiments, this histogram may be used to determine atemporary PTV dose threshold so that only the scores of beamlets withthe highest 50% of PTV doses are considered. This threshold value may bedenoted D_(threshold) ^(PTV) and is the dose of the integrated histogrambin of

with height 0.5.

Next, a BEV score histogram, for determining a score threshold, may becalculated using the BEV score bundle section

, restricted to the indices {n|

_(n) _(v) _(,n) _(x) _(,n) _(y) >D_(threshold) ^(PTV)}. These histogramsmay be denoted as

and

. The ratio of vertical height to horizontal height of a point on theintegrated histogram curve (e.g., the ratio of sides of the rectangleformed by the axes and the dotted lines 730 shown in FIG. 7 ) is:

$\begin{matrix}{{{{ratio}(S)} = {{B^{\mathcal{S}}(S)} \cdot \frac{S_{\max} - S_{\min}}{S - S_{\min}}}},} & (4)\end{matrix}$

where

(S) denotes the integrated histogram height at score value S (i.e. wheren_(bin) is the corresponding bin index). This ratio varies monotonicallyfrom +∞ for S=S_(min) to 0 for S=S_(max) where B^(S)=0. The scorethreshold S_(threshold) may be defined to be the value such thatratio(S_(threshold))=0.7.

It should be understood that the score threshold determination methoddescribed above is only an example. Other determination methods may beused according to other embodiments.

FIGS. 8A and 8B show exemplary images of BEV scores for two adjacentdelivery coordinate vertices for a chest wall treatment. Darker shadesrepresent lower score values and brighter shades represent higher scorevalues. Gold colors represent beamlets that pass the region selectioncriterion, and thus are considered as BEV regions. In this example, fromthe BEV perspective, there are two disconnected BEV regions 810 a and820 a in the BEV shown in FIG. 8A, and two disconnected BEV regions 810b and 820 b in the BEV shown in FIG. 8B.

3. BEV Region Connectivity Manifold

Proximity or overlap of regions at neighboring vertices in the deliverycoordinate space may be examined in order to form a complete BEV regionconnectivity manifold. The BEV region connectivity manifold containsinformation on how candidate target regions in the BEV change, appear,split, and vanish as one moves in the delivery coordinate space in alldirections. For instance, in the example shown in FIGS. 8A and 8B,overlaying the two images of the two adjacent delivery coordinatevertices on each other, the regions 810 a and 810 b on the left mayoverlap with each other, and the regions 820 a and 820 b on the rightmay overlap with each other. Thus, it may be inferred that the regions810 a and 810 b are connected, and the regions 820 a and 820 b areconnected when moving along this edge in delivery coordinate space. Theset of all regions and all connections along delivery coordinate spaceedges form the BEV region connectivity manifold. Moreover, each BEVregion connectivity manifold maintains, for each region it contains, amapping between the region and the subvolumes (that is, voxels in thecomputer model) of the planning target volumes (PTVs), which aresubjected to radiation if the patient body is irradiated from the vertexin which the region resides through an aperture equal to the shape ofthe region.

The information contained in the BEV regions and the BEV regionconnectivity manifold may be used to generate optimized trajectories orfield geometries in a radiation treatment. For example, as discussedbelow and in Christopher Barry Locke and Karl Kenneth Bush, TrajectoryOptimization In Radiotherapy Using Sectioning (TORUS). Medical Physics,2017, methods of trajectory optimization may use a BEV regionconnectivity manifold as a scaffold to guide an optimizer, which maymake the search space small enough to apply graph search techniques withefficient computation times.

IV. Dose Calculations

As discussed above, to determine BEV regions and BEV region connectivitymanifold, it may be necessary to evaluate dose response

_(n) _(v) _(,n) _(x) _(,n) _(y) _(,n) _(ROI) , as expressed in Equation(1), for each point in the 4D space coordinates (n_(v), n_(x), n_(y),n_(ROI)). The dose response may be evaluated by solving the Boltzmanntransport equations. Two approaches of solving the Boltzmann transportequations: the forward transport approach and the adjoint transportapproach, are discussed below. The adjoint transport approach may bemore efficient computationally, as discussed below.

A. Forward Transport Solutions

Evaluating dose response in a forward transport manner may involvesolving the following forward transport equations (see John McGhee, et.al., “AcurosXB Technical Manual,” Varian Medical Systems, (2017)):

{circumflex over (Ω)}·∇Ψ^(γ)+σ_(t) ^(γ)Ψ^(γ) =S ^(γγ)Ψ^(γ) +S ^(γe)Ψ^(γ)+q ^(γ) ,r∈V _(pat),  (5)

{circumflex over (Ω)}·∇Ψ^(e)+σ_(t) ^(e)Ψ^(e) =S ^(ee)Ψ^(e) +S ^(γe)Ψ^(e),r∈V _(pat),  (6)

where Ψ^(γ) is the angular photon fluence, Ψ^(e) is the angular electronfluence, σ_(t) ^(γ) is the photon total cross section, σ_(t) ^(e) is thetotal electron cross section, and S^(γγ), S^(γe), and S^(ee) are thetruncated spherical harmonics source operators. q^(γ) is a point sourcelocated at the position r_(v) ^(s) corresponding to vertex v. Thedependence of Ψ on position r, space angle {circumflex over (Ω)}, andenergy E has been suppressed for clarity. The angular representation forq^(γ) may be defined by a plane of N_(x)×N_(y) pixels centered at theisocenter, where the angle is defined by a vector from the point sourceposition r_(v) ^(s) to the center of each pixel (x_(v,i), y_(v,i)). Forthe purposes of this disclosure, a fluence flow (or intensity) value ofunity may be applied through each pixel. Once Equations (5) and (6) areconsecutively solved with Equation (5) being solved first, then the doseto a region of interest (ROI) may be evaluated as,

D _(ROI)=<σ_(ED) ^(e),Ψ^(e)>≡∫₀ ^(∞)∫_(4π)σ_(ED)Ψ^(e) d{circumflex over(Ω)}dE,r∈V _(ROI),  (7)

Equations (5) and (6) may be discretized by discontinuous finite element(DFEM) in space, multigroup in energy, and discrete ordinates in angle,in conjunction with the first scattered distributed source (FSDS)approach for the point source of photons (see John McGhee, et. al.,“AcurosXB Technical Manual,” Varian Medical Systems, (2017)).

The task of evaluating the dose response using the forward transportapproach can be computationally expensive. As an example, in the contextof a C-arm linear accelerator, assuming that the DCS has resolution of15 degrees in gantry angle and 7.5 degrees in couch angle, the number ofvertices may equal to N_(v)=(360 deg)/(15 deg)×(180 deg)/(7.5 deg)≈580vertices (i.e., candidate directions of incidence). Assuming anisocenter plane size of 40 cm×40 cm and a pixel size of 2.5 mm×2.5 mm,the number of pixels in each isocenter plane may be equal toN=N_(x)×N_(y)=(400/2.5)²=25600. Therefore, in order to characterize thedose for each beamlet separately, the total number of dose calculationsneeded may amount to about N×N_(v)=25600×580≈15×10⁶. Thus, the requiredamount of computation can be very time consuming. For example, thecomputation may take several days, which may render it unsuitable forreal-time treatment optimization in clinical settings.

B. Adjoint Transport Approach

According to some embodiments, the adjoint formulation of the Boltzmanntransport equation may be used for efficient, parallel GPU-compatibleevaluation of the dose matrices. Denoting {tilde over (Ψ)}^(γ) as theadjoint gamma fluence and {tilde over (Ψ)}^(e) as the adjoint electronfluence, the adjoint transport equations may be expressed as,

−{circumflex over (Ω)}·∇{tilde over (Ψ)}^(e)+σ_(t) ^(e){tilde over(Ψ)}^(e) ={tilde over (S)} ^(ee){tilde over (Ψ)}^(e) +{tilde over (S)}^(eγ){tilde over (Ψ)}^(e) +{tilde over (q)} ^(e) ,r∈V _(pat),  (8)

−{circumflex over (Ω)}·∇{tilde over (Ψ)}^(γ)+σ_(t) ^(γ){tilde over(Ψ)}^(γ) ={tilde over (S)} ^(γγ){tilde over (Ψ)}^(γ) +{tilde over (S)}^(γe){tilde over (Ψ)}^(e) ,r∈V _(pat),  (9)

where,

{tilde over (q)} ^(e)=σ_(ED) ,r∈V _(ROI).  (10)

Equations (8) and (9) may be solved consecutively. To obtain the adjointphoton fluence {tilde over (Ψ)}^(γ) at the forward source position r_(v)^(s), the “last collided fluence approach” may be used. Let us assumethat the aperture around the pixel is small enough that only onecharacteristic ray needs to be considered. This characteristic ray maystart at r_(v) ^(s), go through the center of the source fluence pixel(n_(x), n_(y))_(v), and then into the patient (along this definedcharacteristic line), and through the patient and continuing until itexits the patient. A path, s, may be defined to be this characteristicline starting at s_(l)=r_(v) ^(s) and ending at s_(u), where it exitsthe patient. The “last collided fluence approach” uses integraltransport theory to reconstruct the adjoint photon fluence {tilde over(Ψ)}^(γ) at any arbitrary position in space r_(v) ^(s), (see E. E. Lewisand W. F. Miller, Computational Methods of Neutron Transport, AmericanNuclear Society (1993)). The adjoint photon fluence {tilde over (Ψ)}^(γ)at r_(v) ^(s) may be obtained as,

$\begin{matrix}{{{\overset{\sim}{\Psi}}^{\gamma}\left( r_{v}^{s} \right)} = {- {\int_{S_{u}}^{S_{l}}{\left( {{{\overset{˜}{S}}^{\gamma\gamma}{\overset{\sim}{\Psi}}^{\gamma}} + {{\overset{˜}{S}}^{e\gamma}{\overset{\sim}{\Psi}}^{\gamma}}} \right)e^{{- \sigma_{t}}s}{{ds}.}}}}} & (11)\end{matrix}$

In the adjoint transport approach, the dose to a region of interest(ROI), from a unit fluence of photons through one pixel (n_(x),n_(y))_(v) on the point source plane and where the photons originatefrom a single position of source vertex v, may be evaluated as,

D _(ROI) =<q ^(γ),{tilde over (Ψ)}^(γ)>≡∫₀ ^(∞)∫_(4π) q ^(γ){tilde over(Ψ)}^(γ) d{circumflex over (Ω)}dE,r=r _(v) ^(s).  (12)

Note that, for forward transport calculations, D_(ROI) is evaluated overthe volume of the region of interest. For adjoint transportcalculations, D_(ROI) is evaluated at the position of the forward pointsource, r_(v) ^(s).

The adjoint transport approach may afford more efficient doseevaluations as compared to the forward transport approach. As discussedabove, in the forward transport approach, in order to evaluate the dosefrom a unit fluence through each pixel location for each vertex v of theset of all vertices {v_(j)}_(j=1) ^(N) ^(v) , N_(x)×N_(y)×N_(v)×N_(ROI)number of forward computations may be needed with full patient meshtransport solves of Equations (5) and (6), which can be very expensivecomputationally. In contrast, in the adjoint approach, only N_(ROI)adjoint transport solves of Equations (8) and (9) may be needed, alongwith N_(x)×N_(y)×N_(v) ray traces to solve Equation (11). Since Equation(11) can be evaluated independently for a set of vertices in a DCS, thecomputations can be carried out in parallel. Therefore, the number ofexpensive transport solves may be reduced by several orders of magnitudein the adjoint transport approach as compared to the forward transportapproach. A ray trace solve may take only a small fraction of the CPUtime of a transport solve. Thus, the reduction in CPU time may be on theorder of several orders of magnitude as compared to that using theforward transport approach. As a result, the dose evaluations may beperformed in a relatively short time (e.g., in a few minutes to a fewtens of minutes), making it suitable for real-time treatmentoptimization in clinical settings.

C. Dose Evaluations for Multiple Energy Modes

For external beam radiation treatment, various voltages ranging from 4MV to 25 MV may be used for the radiation source (e.g., the voltagesupplied to a linear accelerator in a treatment head). For each voltage,there is an energy spectrum of photons (e.g., X-rays) associated withit. The maximum energy for each voltage may correspond approximately tothe voltage value. For example, for a 6 MV voltage, the maximum energyof the photons may be about 6 MeV.

For existing dose evaluation algorithms using the forward transportapproach, each voltage may require a separate computation of a completetransport. Since transport solves for a single energy in the forwardtransport approach are already very expensive computationally, the taskmay be even more formidable if multiple energies need to be considered.

According to some embodiments, multiple candidate energy modes may beconsidered in a trajectory optimization or beam angle optimization usingthe adjoint transport approach to evaluate dose metrics. In the adjointtransport approach, if Equations (8) and (9) are solved for the maximumenergy possible (e.g., 25 MeV), then any voltage that is less than orequal to 25 MV may be computed without a separate transport solve. Thisis because, if {tilde over (Ψ)}^(γ) is known for all energies up to themaximum energy, any forward beam source spectra, q^(γ), may be used inEquation (12) to obtain D_(ROI). Thus, D_(ROI) can be evaluated for anyenergy less than or equal to the maximum energy.

V. Trajectory Optimization Using Adjoint Transport for Dose

According to some embodiments, trajectory optimization uses the adjointtransport approach to efficiently evaluate dose response of each regionof interest (ROI), for a given direction of incidence (e.g., at vertexn_(v)) in the DCS, to each beamlet at a respective pixel (n_(x), n_(y))of a corresponding BEV plane, as expressed in Equation (12). In someembodiments, a BEV score for each pixel,

_(n) _(v) _(,n) _(x) _(,n) _(y) , may be evaluated as a weighted linearcombination of the doses to the N_(ROI) regions of interest according toEquation (2). BEV regions and a BEV region connectivity manifold may bedetermined based on the BEV scores, e.g., as described above in SectionIII. The information contained in BEV regions and a BEV regionconnectivity manifold may be used to generate optimized trajectories oroptimized field geometries in a radiation treatment, as discussed inmore detail below).

Christopher Barry Locke and Karl Kenneth Bush, Trajectory OptimizationIn Radiotherapy Using Sectioning, Medical Physics, 2017, discussestrajectory optimization in radiotherapy using sectioning (referred toherein as TORUS). The TORUS methods use the BEV regions and BEV regionconnectivity manifold as a guide to generate heuristically optimalradiotherapy trajectories automatically for efficient delivery of highquality VMAT treatment plans. TORUS uses an optimization graph on top ofa delivery coordinate space to generate optimal treatment trajectoriesusing a dual-metric optimization. Nodes in the optimization graph mayrepresent individual control points, and trajectories may be defined aspaths that minimize a min-distance metric, while a max-distance metricmay act as a measure of goodness to select optimal trajectories.

A. PTV Angular Flux

One of the concepts used in the TORUS methods is PTV angular flux, whichrelates to novelty of three-dimensional (3D) direction vectors ofincident beamlets for a given point in a PTV. Inverse dose optimizationmay perform better with more angles from which radiation beams enter thepatient. The reason for this may be twofold. First, by entering thepatient from many directions, the ratio of overlapping dose within thePTV to surrounding OAR may be greater, resulting in steeper dosegradients outside the PTV. Second, each beamlet from each direction mayprovide a different 3D dose contribution to the patient. Therefore,increasing the number of such unique beamlets may give the optimizermore “basis vectors” to work with when sculpting optimal dose profilesaround critical structures.

Note that just entering the patient from many directions may not besufficient to ensure optimal plan quality. It may be that in some cases,portions of the PTV are only able to be exposed from a small number ofdirections when protecting nearby OAR, even though the number of beamsentering the patient is high. This can result in either under coverageof small regions of the PTV, non-conformal regions (dose streaks), orunsatisfactory dose compensation. To encourage maximal coverage andconformality, it may be desirable to have each elemental volume of thePTV be individually targeted from many different directions.

According to some embodiments, an angular flux of a given point in a PTVmay be evaluated by computing the 3D direction vectors of incidentbeamlets and binning them in angular bins. FIG. 9 illustrates anexemplary quadrilateralized spherical cube that may be used forevaluating the angular flux. As illustrated, each cube face may bedivided into 4×4 squares. Thus, there are 6×4×4=96 squares. Each squarecorresponds to a single angular bin. This may provide a bin size on theorder of 20 degrees, which may correspond to the same order of distancebetween vertices in a delivery coordinate space (DCS). In a generalcase, each cube face may be divided into 2^(2n) squares, where n is apositive integer. In the example illustrated in FIG. 9 , n=2. Thisbinning method may result in bins of unequal solid angles, but thedifferences may be relatively small (at most 19%). Additionally,randomly orienting each spherical cube across a set of sampling pointsmay cancel out such discrepancies. In other embodiments, the angularflux may be evaluated by computing the 3D direction vectors of incidentbeamlets through a closed surface of a shape other than thequadrilateralized spherical cube. For example, the closed surface can bea spherical surface or a cube surface.

A set of sampling points may be distributed within a PTV. The angularflux at each sampling point may be evaluated and optimized duringtrajectory optimization. According to some embodiments, indexing eachbin from 0 to 95, the angular flux of a given PTV point may be stored asa 12 byte bitset, enabling fast bitwise calculations to be performed. Ifthere are n PTV points considered, then the angular flux state

is represented as a vector of length n of bitsets each with 96 bits.

B. Dual-Distance Metric

According to some embodiments, the information stored in the BEV scorebundle sections and the BEV region connectivity manifold may be used togenerate treatment plans for radiotherapy. The BEV region connectivitymanifold may serve as a scaffold to guide the optimizer, which may makea search space small enough to apply graph search techniques with fastcomputation times.

In general, optimal trajectories in a VMAT-like treatment may be onesthat hit the PTV as much as possible, avoid or minimize healthy tissuedoses, and enter the PTV from many different directions, and can becompleted in a relatively short delivery time. Some of the goals mayconflict with each other. For example, a treatment plan that treats eachPVT element from every direction with a high degree of MLC modulationmay produce a nearly ideal dose distribution, but may also take anexcessive amount of time to deliver. Therefore, it may be desirable tofind a trajectory that covers many directions around the patient in anefficient way. According to some embodiments, an optimization method mayseek to maximize trajectory length, while using relatively “straight”trajectories for delivery efficiently. The “straightness” of atrajectory may be understood in terms of a geodesic line. A geodesic isthe shortest path between two given points in a curved space. A geodesicmay be calculated by finding a line that minimizes a distance functionbetween the two points.

According to some embodiments, based on the BEV score bundle sectionsand the BEV region connectivity manifold, trajectories that have controlpoints traversing through regions of high region scores R may bepreferred (e.g., as defined in Equation (3)). To encourage geodesics totraverse such control points, a distance function may be defined suchthat smaller non-negative distances are preferred. On the other hand, topick long trajectories, another distance function may be defined suchthat larger distances are preferred.

To overcome this inherent conflict between trying to find “short” aswell as “long” trajectories, a dual-distance metric approach is used inan optimization according to some embodiments. The dual-distance metricincludes two distance functions that play different roles in theoptimization. The two distance functions may be referred to as amin-distance function D^(min) (where smaller values are preferred;minimization defines geodesic paths in a graph), and a max-distancefunction D^(max) (where larger values are preferred; to be maximizedthrough selection of trajectory).

C. Stateful Graph Optimization

According to some embodiments, a symmetric directed graph may be usedfor trajectory optimization using the min-distance function D^(min) andthe max-distance function D^(max). Since the distance functions maydepend not just on the edges of the graph, but also on the history of atrajectory up to that point, a statefulness may be introduced to thegraph.

The graph may include a set of nodes, and a set of edges that connectthe nodes. To account for the history in the graph, a point may bedefined to be a node and state pair P=(N,

). The state may be the PTV angular flux state described above. Given apoint P₁=(N₁,

₁) and an edge E=(N₁, N₂), the successor point is P₂=(N₂,

₂), where the new state

₂ may be given by a successor state function σ(

₁, E), which describes how the state changes moving along the edge E.The min-distance and max-distance between these two points may bedenoted as D^(min/max)(P₁, P₂).

A path may be an ordered sequence of points, P(P₁, . . . , P_(n)) suchthat ∃E=(N_(i), N_(i+1)) and

_(i+1)=σ(

_(i), (N_(i), N_(i+1))) for all 1≤i<n. The min-distance and max-distanceof this path may be D^(min/max)(P)=Σ_(i=1) ^(n−1)D^(min)/max(P_(i),P_(i+1)) Given an initial state

and two nodes N and N′, let the set of all possible paths between them,P(P₁, . . . , P_(n)) such that P₁=(N,

) and P_(n)=(N′,

′), where

′ is unspecified, be denoted

(N, N′|

). A set of trajectories may be defined to be the set of paths with aminimal min-distance,

(N,N′|

)=arg

D ^(min)(P)  (13)

A set of optimal trajectories between these two nodes may then bedefined to be the trajectories with a maximum max-distance,

^(optimal)(N,N′|

)=arg

D ^(max)(P).  (14)

Using the definitions provided by Equations (5) and (6), themin-distance and max-distance between two arbitrary nodes isD^(min/max)(N,N′|

)=D^(min/max)(P) for P∈

^(optimal)(N,N′|

).

Given an initial state

, the goal now may be to define a globally optimal trajectory. Ifattention is restricted to start and end nodes in some set

, then a set of optimal trajectories ending on the set

may be defined to be:

(

)=

^(optimal)(N,N′|

)

where

=arg

D _(max)(N,N′|

).  (15)

The optimization problem may be to find an element of

(

) for some initial state

and the set of start and end points

.

D. Graph Definitions

According to some embodiments, a control point in the graph may beuniquely determined by three integers, vertex v (i.e., a point on thedelivery coordinate space or DCS), collimator index c (which determinescollimator angle out of a discrete set of possibilities), and regionbitfield b. The region bitfield b is a list of boolean flags thatdetermine which subset of regions to select for a given vertex. Thesethree integers may define a node in the graph as N=(v, c, b). StartingMLC leaf positions may be determined by fitting to this subset ofregions of the BEV.

According to some embodiments, the BEV region connectivity manifold mayinclude multiple mutually disjoint connected components, each of whichforms a single connected part of the total search graph through thefollowing definition of edges. Given two nodes N₁=(v₁, c₁, b₁) andN₂=(v₂, c₂, b₂), there may be an edge E connecting these two nodes ifthe following are satisfied:

-   -   There is an edge connecting v₁ and v₂ in the delivery coordinate        space;    -   Δt_(collimator)≤Δt_(directional), where Δt_(collimator) is the        time for the collimator to move        Δθ_(collimator)=θ_(collimator)(c₂)−θ_(collimator)(c₁), and        Δt_(directional) is the time to move in gantry and couch space        from vertex v₁ to v₂.

The set of boundary nodes

are the potential start and end nodes in the graph optimization. Todefine this set of boundary nodes, it may be necessary to first defineboundary vertices in the delivery coordinate space for one-dimensional(1D) and two-dimensional (2D) spaces (which may be generalized to higherdimensions). A 1D space may be made of only vertices and edges, andboundary vertices are those which touch at most a single edge.Similarly, in a 2D space of vertices, edges and faces, boundary verticesare those which belong to an edge that only touches a single face. Withthis definition, a region

=(n_(v), n_(subindex)) may be a boundary region if one of the followingconditions is satisfied:

-   -   (1) Vertex v of index n_(v) is a boundary vertex;    -   (2) There exists an edge in the delivery coordinate space        touching vertex v of index n_(v) such that there is no region        edge along this delivery coordinate space edge emanating from        region r.

FIGS. 10A-10B illustrate these conditions. FIG. 10A may represent a 1Ddelivery coordinate space, with boundary vertices 1010 and 1020represented as filled circles. FIG. 10B may represent a 2D deliverycoordinate space. The BEV region connectivity manifold is represented byovals (with each oval representing a region) in five BEV planes 1030a-1030 e, and the connections between them. The region 1040 in theleftmost plane 1030 a and the region 1050 on the rightmost plane 1030 eare boundary regions that satisfy condition (1) above, while the region1060 in the plane 1030 d has no connection to a region on the right,thus is a boundary region that satisfy condition (2) above. With thesedefinitions, a node is a boundary node if all the regions denoted by itsregion bitfield are boundary regions.

The state information used in the graph optimization is PTV angular flux

. Using a bitset definition of PTV angular flux as described above, thesuccessor function σ, which defines how the state changes, may bedefined by the bitwise OR operator, σ(

₁, (N₁, N₂))=

₁|

(N₂), where

(N₂) is the contribution to the angular flux state from the regionsdenoted by region bitfield b₂ at vertex v₂.

Given two points P_(i)=(N_(i),

_(i)) where N_(i)=(v_(i), c_(i), b_(i)) for i=1, 2 with edge Econnecting the nodes, min/max-distance functions may be defined as,

$\begin{matrix}{{{D^{\min}\left( {P_{1},P_{2}} \right)} = {\Delta{\theta(E)}\left( {\frac{3}{S\left( {P_{1},P_{2}} \right)} + {t^{mlc}(E)} + {p^{{eff}.}(E)}} \right)}},} & (16)\end{matrix}$ $\begin{matrix}{{{D^{\max}\left( {P_{1},P_{2}} \right)} = {\Delta{\theta(E)}{S\left( {P_{1},P_{2}} \right)}}},} & (17)\end{matrix}$

where Δθ is the physical angular distance traveled by the treatmenthead. The min-distance function is set to +∞ if the resulting MLCconfiguration violates machine limitations. Note that the score termS(P₁, P₂) appears as the inverse of one another in each equation,reflecting the fact that, roughly speaking, max-distance=“goodness” andmin-distance=1/“goodness”. The definitions and meanings behind each termare as follows:

-   -   S(P₁, P₂)=S^(score)(b₂)+2S^(angular)(        ₁,        ₂)

${{S^{score}(b)} = {\frac{{\mathcal{R}_{avg}(b)} + 1}{2} \cdot {\max\left( {0,{{A(b)} - {A_{penalty}(b)}}} \right)}}},$

-   -    where        _(avg)(b) is the average region score of the regions of bitfield        b per unit area, A(b) is the combined area of regions of        bitfield b, and A_(penalty)(b) is the total non-region area        exposed by the fitted MLCs. This term penalizes poor MLC target        fitting and encourages high scoring regions.    -   S^(angular)(        ₁,        ₂) is a term that encourages regions that provide novel angles        to the existing angular flux state. If        _(avg) is the average contribution (number of bits) to a blank        angular flux state from each region, then S^(angular)(        ₁,        ₂) is the number of bits in (        ₂&˜        ₁) or 1 (whichever is larger), normalized by dividing by        _(avg).    -   t^(mlc) is the time for the MLCs to move between control points.        This penalizes collimator angles that result in excessive MLC        motion.    -   p^(eff.)=(1−Δθ/Δθ_(max)) is a factor that penalizes edges where        the trajectory is almost stationary, so treatment time is not        wasted in such locations.

The score term S(P₁, P₂) may be the main force driving the optimizer tofind trajectories that target the PTV from good directions and that givecontributions from different directions while avoiding poor MLC targetfitting.

E. Graph Optimization Solution

The stateful dual-distance metric graph optimization defined above maybe solved using the Dijkstra algorithm. Converting the list of pointsfrom the resulting optimal trajectory to control points may provide thedesired radiotherapy trajectory. Running Dijkstra's algorithm from astart node without stopping at any particular end node may result in atree structure of points, effectively completing a 1-to-N search fromthe given starting node to all other nodes with the same computationalcomplexity as the usual 1-to-1 search between two nodes. By selectingthe trajectory with the largest max-distance, the restrictedoptimization problem of finding the optimal trajectory from a given nodemay be solved efficiently.

In general, to find the globally optimal trajectory with the largestmax-distance may require repeating this computation from every possiblenode (N-to-N search). An approximately optimal trajectory can be foundby picking an arbitrary start node, and running this algorithmrepeatedly with the same initial state, using the end node of theprevious run as the start node of the successive run. In someembodiments, this process is repeated twice; thus trajectoryoptimization may be carried out with the same computational complexityas the underlying Dijkstra algorithm.

The path optimization may be carried out for each connected component ofthe graph. The trajectory with the largest max-distance across allpossibilities may be selected in the end. The presence of the angularflux state in the distance function may help ensure that the selectedtrajectory will also be one that tends to provide novel directions fromwhich to treat the PTV.

F. Method of Trajectory Optimization for Radiotherapy Treatment UsingSectioning

FIG. 11 shows a flowchart illustrating a method 1100 of trajectoryoptimization for radiotherapy treatment using sectioning according tosome embodiments.

At 1102, a patient model is provided. The patient model includes one ormore regions of interest (ROIs) for the radiotherapy treatment.

At 1104, a delivery coordinate space (DCS) is defined. The DCS has a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane.

At 1106, for each respective ROI of the one or more ROIs, an adjointtransport is solved to obtain an adjoint solution field from therespective ROI. For each respective candidate vertex in the DCS, and foreach respective pixel of the respective BEV plane defined by therespective candidate vertex, an adjoint photon fluence originating froma respective beamlet incident from the respective candidate vertex andpassing through the respective pixel is evaluated by performing raytracing of the adjoint solution field; and a respective dose of therespective ROI from the respective beamlet is evaluated using theadjoint photon fluence.

At 1108, for each respective candidate vertex in the DCS, and for eachrespective pixel of the respective BEV plane defined by the respectivecandidate vertex, a respective BEV score of the respective pixel isevaluated using the doses of the one or more ROIs evaluated for therespective beamlet incident from the respective candidate vertex andpassing through the respective pixel; and one or more BEV regions in therespective BEV plane are determined based on the BEV scores of thepixels of the BEV planes corresponding to the set of candidate verticesin the DCS.

At 1110, a BEV region connectivity manifold is determined based on theBEV regions of the BEV planes of the set of candidate vertices in theDCS. The BEV region connectivity manifold represents connections betweencontiguous BEV regions between adjacent vertices.

At 1112, one or more optimal treatment trajectories are determined basedon the BEV region connectivity manifold.

It should be appreciated that the specific steps illustrated in FIG. 11provide a particular method of trajectory optimization according to someembodiments. Other sequences of steps may also be performed according toalternative embodiments. For example, alternative embodiments mayperform the steps in a different order. Moreover, the individual stepsmay include multiple sub-steps that may be performed in varioussequences as appropriate to the individual step. Furthermore, additionalsteps may be added and some steps may be removed depending on theparticular applications. One of ordinary skill in the art wouldrecognize many variations, modifications, and alternatives.

VI. Beam Angle Optimization Using Adjoint Transport for Dose

The TORUS methods use the BEV regions and the BEV region connectivitymanifold as a guide to generate heuristically optimal radiotherapytrajectories automatically for efficient delivery of high quality VMATtreatment plans. According to some embodiments, the TORUS methods aremodified to generate heuristically optimal IMRT fields. As discussedabove, the TORUS methods use an optimization graph on top of a deliverycoordinate space to generate optimal treatment trajectories through theuse of a dual-metric optimization. Nodes in the optimization graphrepresent individual control points, and trajectories are defined to bepaths that minimize a min-distance metric, while a max-distance metricacts as a measure of goodness to select optimal trajectories. For anIMRT treatment plan, treatment fields may be represented in thisframework as a set of nodes. A beam angle optimization may involvefinding a set of k nodes that have optimal max-distance.

A. Field Geometry Optimization

According to some embodiments, to find a set of k beams, a BAO graphbased on the TORUS graph concept may be built based on the followingnodes and edges:

Nodes

Nodes are a set

of N tuples of the form (v, c, b), where v is a vertex, c is thecollimator angle index, and b is a binary mask of included regions. Eachtuplet may determine a single field. The vertex v may correspond to alocation in a discretized delivery coordinate space (DCS). For example,in a C-arm linear accelerator, the vertex v may correspond to anisocenter, a gantry angle θ_(gantry), and a couch angle θ_(couch). Thecollimator angle index c may correspond to a collimator angle θ_(c) outof a set of discrete possible collimator angles. The region mask b maycorrespond a set of contiguous target regions the MLC leaves may expose.

Edges

Edges are connections between neighbor nodes that have vertex-vertexconnectivity in the underlying delivery coordinate space, and MLCconnectivity between the respective MLC leaf sequences. In the case ofbeam angle optimization, the MLC connectivity constraint may not beimportant, but it may nevertheless act to reduce the number of edges inthe graph. Thus, it may be computationally useful to keep the MLCconnectivity constraint.

In the TORUS methods, a min-distance and a max-distance may be definedalong edges in the TORUS graph. In the case of static fields, there isno gantry motion while the treatment beam is on. Therefore, only thevertices themselves need to be considered, and there is no need for amin-distance function. The score (max-distance) may be the metric tooptimize. The score S may be defined on subsets of nodes as,

S(

)=max−distance corresponding to the set of nodes (beams)

⊂

.  (18)

The optimal set of k beams may be defined to be the subset

⊂

with |

|=k that gives an optimal score S(

), where |

| represents the number of beams in the subset

⊂

.

The score S(

) may be a complex non-local function of the entire subset

. Therefore, it may be a non-trivial problem to solve exactly. However,there may exist an efficient way to find approximate solutions fasterthan trying all _(N)C_(k) combinations. In some embodiments, a beamangle optimization may include the following steps:

-   -   (1) Evaluate the scores for individual beams (|        |=1), and use the squares of these scores as sampling        probabilities.    -   (2) Randomly sample k beams using the sampling probabilities.    -   (3) Apply local gradient descent to this set        to find a local minimum.    -   (4) Repeat steps (1)-(3) until there is no improvement in score        found for a predetermined number of successive trials (e.g., 20        successive trials).

This procedure may be referred to as a grandient descent method. Duringeach iteration of the gradient descent procedure, all possible “edges”may be considered looking for improvement, where an “edge” is motion ofthe subset

to a neighbor subset

′. The subset of beams may be updated to move in the direction that mostimproves the score. This procedure may be repeated until no more localimprovements in the score is found. A neighbor to subset

is defined to be another subset

′ that differs in exactly one beam B_(i)→B′_(i), where these two beamseither share an edge in the TORUS graph or have the same position in thedelivery coordinate space (same vertex).

B. Score Function

The score function S(

) may be similar to the max-distance function as in TORUS. The score S(

) of a set of beams

may include two parts: individual beam scoring, and overall PTV angularflux novelty. The presence of the global PTV angular flux metric is whatmakes the score function non-local, whereby it is not simply a functionof the scores of the individual beams. The score S(

) may be written as the sum of two parts,

S(

)=S _(local)(

)+3S _(flux)(

),  (19)

where the local part S_(local)(

) may be defined as

$\begin{matrix}{{{S_{local}(B)} = {{\sum}_{B \in B}{S_{beam}(B)}}},} & (20)\end{matrix}$ $\begin{matrix}{{{S_{beam}(B)} = {{w(B)} \cdot \left( {\frac{s(B)}{A_{avg}} + \frac{0.2}{\max\left( {0.05,{o(B)}} \right)}} \right)}},} & (21)\end{matrix}$ $\begin{matrix}{{{w(B)} = {\max\left( {{0\text{.1}},{\min\left( {{1.0},{{1.1} - {c^{2}(B)}}} \right)}} \right)}},} & (22)\end{matrix}$ $\begin{matrix}{{{c(B)} = {{MLC}{contention}{severity}}},} & (23)\end{matrix}$ $\begin{matrix}{{{s(B)} = {{integrated}{score}{of}{regions}}},} & (24)\end{matrix}$ $\begin{matrix}{{o(B)} = {\frac{2}{1 + {{CrossSectio}{n^{2}(B)}}} + {{x_{rms}(B)}/2.}}} & (25)\end{matrix}$

A_(avg) may be defined to be the average area of regions. Thus, the term

$\frac{s(B)}{A_{avg}}$

may be considered as the average region score. In some embodiments, theMLC contention severity may relate to A_(penalty) defined above as thetotal non-region area exposed by the fitted MLCs. The cross section maybe defined to be the y extent of the open MLC leaves normalized to theaverage diameter of regions as circles. x_(rms)(B) may be defined as theroot-mean-square x extent of open leaves. The flux part S_(flux)(

) may be defined as the PTV angular flux novelty of the angular fluxstate generated by the selected regions, normalized by dividing by theaverage angular flux novelty of the individual regions. The separationof the score into local and non-local parts may allow some codeoptimizations to be done by pre-computing the contributions ofindividual beams to the local part.

C. Beam Angle Optimization Including Consideration of Beam-Off Time

Consider an IMRT treatment as an external-beam radiotherapy treatment,where the dose is delivered from k static beam locations r_(i)=(v_(i),c_(i)), i=1, . . . , k, where the indices v_(i) and c_(i) correspond tovertex and collimator angle index, respectively. In a C-arm linearaccelerator treatment system, a vertex may include an isocenter, agantry angle, and a couch angle. In other types of external-beamradiation treatment systems, a vertex may include other treatment axesvariables. The fluence delivered from each r_(i) may be determined byway of a fluence-optimization scheme that distributes the fluenceoptimally between the r_(i) such that clinical optimization objectivesare fulfilled.

A dosimetrically optimal IMRT treatment plan may be found within a givendelivery coordinate space by running fluence optimization and leafsequencing for each permissible combination of {r_(i)}_(i=1) ^(k), andby picking the dosimetrically optimal one. However, this process may beimpractical due to the large number of possible combinations. Moreover,it may be desirable that the patient and the patient's internal organsremain stationary during the administration of radiation for thedelivered dose to match the planned dose. The longer the treatmenttakes, the more likely that the patient or the patient's internal organsmay move during the treatment, and hence a higher probability of notdelivering the intended dose to the target volumes.

The overall treatment time may be prolonged by the beam-off transitiontimes from beam r_(i) to r_(i+1). The total beam-off transition timefrom beam r_(i) to r_(i+1) may be expressed asΔt_(i,i+1)=max_(j)Δt_(i,i+1) ^(change) ^(j) , where the changes mayinclude, but are limited to the following:

-   -   a time of Δt_(i,i+1) ^(move) to move machine axes from positions        at vertex v_(i) to those at v_(i+1);    -   a time of Δt_(i,i+1) ^(collimator) to rotate the collimator from        c_(i) to c_(i+1);    -   a time of Δt_(i,i+1) ^(imaging) of acquiring an image or        multiple images for guidance of treatment between beams i and        i+1.

A minimization of the aggregate beam-off time Δt=_(i=1) ^(k−1)Δt_(i,i+1) may be an important part of the beam-selection problem. Insome embodiments, a BAO algorithm may include constraints on pickingbeams that can be delivered in a time-efficient radiation order. In someother embodiments, the radiation order of the beams may be determined asa post-processing step.

D. Method of Beam Angle Optimization Using Sectioning

FIG. 12 shows a flowchart illustrating a method 1200 of beam angleoptimization using sectioning for an IMRT radiotherapy treatmentaccording to some embodiments.

At 1202, a patient model is provided. The patient model includes one ormore regions of interest (ROIs) for the IMRT radiotherapy treatment.

At 1204, a delivery coordinate space (DCS) is defined. The DCS has a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane.

At 1206, for each respective ROI of the one or more ROIs, an adjointtransport is solved to obtain an adjoint solution field from therespective ROI; and for each respective candidate vertex in the DCS, andfor each respective pixel of the respective BEV plane defined by therespective candidate vertex: an adjoint photon fluence originating froma respective beamlet incident from the respective candidate vertex andpassing through the respective pixel is evaluated by performing raytracing of the adjoint solution field; and a respective dose of therespective ROI from the respective beamlet is evaluated using theadjoint photon fluence.

At 1208, for each respective candidate vertex in the DCS, and for eachrespective pixel of the respective BEV plane defined by the respectivecandidate vertex: a respective BEV score of the respective pixel isevaluated using the doses of the one or more ROIs evaluated for therespective beamlet incident from the respective candidate vertex andpassing through the respective pixel; and one or more BEV regions in therespective BEV plane are determined based on the BEV scores of thepixels of the BEV planes corresponding to the set of candidate verticesin the DCS.

At 1210, a BEV region connectivity manifold is determined based on theBEV regions of the BEV planes of the set of candidate vertices in theDCS. The BEV region connectivity manifold represents connections betweencontiguous BEV regions between adjacent vertices.

At 1212, a set of IMRT fields is determined based on the BEV regionconnectivity manifold. Each respective IMRT field of the set of IMRTfields defines a beam angle corresponding to a respective vertex in theDCS.

It should be appreciated that the specific steps illustrated in FIG. 12provide a particular method of beam angle optimization according to someembodiments. Other sequences of steps may also be performed according toalternative embodiments. For example, alternative embodiments mayperform the steps in a different order. Moreover, the individual stepsmay include multiple sub-steps that may be performed in varioussequences as appropriate to the individual step. Furthermore, additionalsteps may be added and some steps may be removed depending on theparticular applications. One of ordinary skill in the art wouldrecognize many variations, modifications, and alternatives.

VII. Trajectory Optimization Considering Multiple Candidate Energy Modes

As discussed above, for a given direction of incidence (e.g., a givenvertex n_(v) in the DCS), a dose response of a given region of interestn_(ROI) to a beamlet at a respective pixel (n_(x), n_(y)) of acorresponding BEV plane,

_(n) _(v) _(,n) _(x) _(,n) _(y) _(,n) _(ROI) , may be evaluated usingthe adjoint transport approach according to Equation (12). In someembodiments, the dose response may be evaluated using other methods. Forexample, the dose response may be evaluated using the forward transportapproach according to Equation (7). A BEV score for each pixel,

_(n) _(v) _(,n) _(x) _(,n) _(y) , may be evaluated as a weighted linearcombination of the doses to the N_(ROI) regions of interest according toEquation (2). The BEV score may be a measure of “goodness” of thebeamlet as a candidate for irradiation of the target volumes. For agiven BEV plane, BEV regions may be identified as candidate irradiationapertures according to a region selection criterion based on the BEVscores.

The size of a BEV region may depend on all other candidate directions.The size of a BEV region may also depend on how deep-seated the targetvolumes (PTVs) are, and where organs at risk (OARs) are located relativeto the target volumes along the direction of incidence. In addition, thesize of a BEV region may depend on the energy mode of a radiation beam.The energy mode of a radiation beam may be a combination of the photonenergy and the distribution of photon number density perpendicular tothe direction of irradiation in the divergent frame of reference. Thedependence of a BEV region on energy mode may be due to the differencesof depth-dose curves at different energy modes. For example, for a givenphoton number density, at low photon energies, more energy may bedeposited in superficial tissues than at high photon energies; whereasat high photon energies, more energy may be deposited in deep-seatedtissues.

FIGS. 13A and 13B show exemplary BEV regions for a given direction ofincidence for a first energy mode and a second energy mode,respectively, where the second energy mode has a higher energy than thefirst energy mode, according to some embodiments. Darker shadesrepresent lower score values and brighter shades represent higher scorevalues. Orange colors represent beamlets that pass the region selectioncriterion, and thus are considered as BEV regions as potential openaperture candidates for irradiation of target structures. The bluecolored areas correspond to those pixels through which irradiation maydeposit dose in the target structures. The light green colored areascorrespond to those pixels through which irradiation may not passthrough any target structures.

A region score R for a BEV region, in a given direction of incidence,may be defined as,

R=∫R(n _(x) ,n _(y))dA,  (26)

where R(n_(x), n_(y)), a non-negative real number smaller than or equalto one, is a normalized weighted sum of the aggregate dose deposited toregions of interest from a beamlet corresponding to pixel (n_(x),n_(y)), and the integration is over the area of the BEV region (e.g.,the orange colored region in FIG. 13A or FIG. 13B). The larger theregion score R, the more dose is deposited to the target volume and lessdose is deposited to organs at risk via irradiation through the BEVregion.

As illustrated in FIGS. 13A and 13B, a higher energy mode (FIG. 13B) maylead to a larger and brighter BEV region as compared to a lower energymode (FIG. 13A). A larger area of a BEV region means that the targetvolume may be safely irradiated through a larger aperture. A brighterBEV region means that more aggregate dose may be deposited to the targetvolume. Based on Equation (26), both a larger and brighter BEV regionmay lead to a higher region score R. Therefore, because of the possibledependence of BEV regions on energy, a given direction of incidence maybe preferable for some energy modes, while not for some other energymodes.

According to some embodiments, trajectory optimization and beam angleoptimization in a radiation treatment plan may consider multiplecandidate energy modes. For example, instead of considering a singleenergy mode, a set of candidate energy modes may be considered. BEVregions and BEV region connectivity manifold may be identified andconstructed for each energy mode of the set of candidate energy modes.Trajectory optimization may be performed concurrently with theoptimization of energy modes based on the BEV regions and the BEV regionconnectivity manifold for the set of candidate energy modes.

In the following, two approaches of trajectory optimization consideringmultiple candidate energy modes are described. In a first approach, atrajectory optimization is performed for each energy mode of a set ofcandidate energy modes to obtain a respective candidate set oftrajectories for each energy mode. The energy mode corresponding to thecandidate set of trajectories that gives rise to the optimal value of anobjective function (“score”) may be picked. In some embodiments, theobjective function may be the dual max-min objective function in TORUS,as described above. The first approach may be referred to assingle-energy-mode trajectory optimization over multiple candidateenergy modes, as the set of trajectories selected for the radiationtreatment plan has the same energy mode for all trajectories in the set.

In a second approach, the end result of a trajectory optimization may bea set of trajectories, possibly with different energy modes fordifferent trajectories. The second approach may be referred to asmixed-energy-mode trajectory optimization over multiple candidate energymodes, as some trajectories in the set may be in first energy mode whilesome other trajectories in the set may be in a second energy modedifferent from the first energy mode. The capability to mix differentenergy modes may be beneficial in treatments in which a homogeneous doseis a high-priority clinical goal, or in treatments in which some targetsare superficial while some other targets are deep-seated (e.g., inintracranial multi-metastasis cases or breast multi-metastasis cases).

A. Single-Energy-Mode Trajectory Optimization Over Multiple CandidateEnergy Modes

In some embodiments, a methods of trajectory optimization considering Mcandidate energy modes may generate a set of N trajectories intended forVMAT-like delivery using a single energy mode by running a pathoptimization algorithm sequentially M times, as described in an examplebelow.

According to some embodiments, to find an optimal energy mode among Mcandidate energy modes (indexed as m=1, . . . M) for a set of Ntrajectories (indexed as n=1, . . . N), a trajectory optimization mayinclude the following steps:

-   -   1. Set m=1;    -   (a) Set n=1 and PTV angular flux state F₀ ^(m)=0;    -   (i) Find n^(th) optimal trajectory T_(n) ^(m) and flux state        F_(n) ^(m) against flux state F_(n−1) ^(m) at energy mode m.    -   (ii) If n=N, go to 1(b). Otherwise, set n=n+1 and go to 1(a)(i).    -   (b) If the set of trajectories {T_(n) ^(m)}_(n=1) ^(N) is better        than the current energy mode optimum set        , based on an objective function, set        _(current)={T_(n) ^(m)}_(n=1) ^(N) and F_(current)=F_(N) ^(m);    -   2. If m=M, stop. Otherwise, set m=m+1 and go to 1(a).

In some embodiments, considering VMAT-like trajectories, along which oneor several machine axes move while the treatment beam is on, ifconnected BEV regions are large for a given energy mode in a set of DCSvertices connected by edges, that energy mode may be considered as agood candidate energy mode, and the set of DCS vertices may beconsidered as a good candidate path for a VMAT trajectory using thatenergy mode.

B. First Method of Trajectory Optimization Considering MultipleCandidate Energy Modes

FIG. 14 is a flowchart illustrating a method 1400 of trajectoryoptimization considering multiple candidate energy modes according tosome embodiments. The method 1400 may be referred to assingle-energy-mode trajectory optimization over multiple candidateenergy modes, as discussed above.

At 1402, a patient model is provided. The patient model includes one ormore regions of interest (ROIs) for the radiotherapy treatment.

At 1404, a delivery coordinate space (DCS) is defined. The DCS has a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane.

At 1406, a plurality of candidate energy modes are identified for theradiotherapy treatment.

At 1408, for each respective energy mode of the plurality of candidateenergy modes, for each respective BEV plane of a respective candidatevertex in the DCS, and for each respective ROI of the one or more ROIs:a respective dose of the respective ROI from a respective beamletincident from the respective candidate vertex and passing through arespective pixel of the respective BEV plane is evaluated usingtransport solutions for the respective energy mode; a respective BEVscore of the respective pixel is evaluated using the doses of the one ormore ROIs evaluated for the respective beamlet; one or more BEV regionsin the respective BEV plane for the respective energy mode aredetermined based on the BEV scores of the pixels of the BEV planescorresponding to the set of candidate vertices in the DCS; and arespective BEV region connectivity manifold for the respective energymode is determined based on the BEV regions of the BEV planes of the setof candidate vertices in the DCS. The respective BEV region connectivitymanifold represents connections between contiguous BEV regions betweenadjacent vertices.

At 1410, a plurality of candidate sets of optimal treatment trajectoriesis determined by: for each respective energy mode of the plurality ofcandidate energy modes, determining a respective candidate set ofoptimal treatment trajectories based on the respective BEV regionconnectivity manifold for the respective energy mode.

At 1412, one of the plurality of candidate sets of optimal treatmenttrajectories is selected as a final set of optimal treatmenttrajectories based on an objective function. The final set of optimaltreatment trajectories corresponds an optimal energy mode among theplurality of candidate energy modes.

It should be appreciated that the specific steps illustrated in FIG. 14provide a particular method of trajectory optimization according to someembodiments. Other sequences of steps may also be performed according toalternative embodiments. For example, alternative embodiments mayperform the steps in a different order. Moreover, the individual stepsmay include multiple sub-steps that may be performed in varioussequences as appropriate to the individual step. Furthermore, additionalsteps may be added and some steps may be removed depending on theparticular applications. One of ordinary skill in the art wouldrecognize many variations, modifications, and alternatives.

C. Mixed-Energy-Mode Trajectory Optimization Over Multiple CandidateEnergy Modes

According to some embodiments, to find a set of N optimal trajectoriesamong M candidate energy modes, in which the N trajectories may havedifferent energy modes (referred to herein as mixed energy mode), atrajectory optimization considering multiple candidate energy modes mayinclude the following steps:

-   -   1. Set n=1 and PTV angular flux state F_(current)=0;    -   (a) Set m=1;    -   (i) Find n^(th) optimal trajectory T_(n) ^(m) and flux state        F^(m) against flux state F_(current) at energy mode m.    -   (ii) If T_(n) ^(m) is better than the current optimum        T_(current), based on the objective function, set        T_(current)=T_(n) ^(m) and F_(n)=F^(m);    -   (iii) If m=M, go to 1(b). Otherwise, set m=m+1 and go to        1(a)(i).    -   (b) Set T_(n)=T_(current) and F_(current)=F_(n), where        F_(current) reflects the PTV angular flux state after finding n        trajectories.    -   2. If n=N, stop. Otherwise, set n=n+1 and go to 1(a).

In some embodiments, a constraint of at most one energy mode change canbe added. That is, if n=1, . . . , k trajectories have a constant energymode m₀, but energy mode m₁≠m₀ is found for trajectory n=k+1, then theenergy mode m₁ is kept for the remaining n=k+2, . . . , N trajectories.

In some embodiments, considering VMAT-like trajectories, along which oneor several machine axes move while the treatment beam is on, ifconnected BEV regions are large for a first energy mode in a first setof DCS vertices connected by edges, and connected BEV regions are largefor a second energy mode in a second set of DCS vertices connected byedges, the first set of DCS vertices and the second set of DCS verticesmay be considered as a first candidate trajectory and a second candidatetrajectory, respectively, for a mixed-energy mode VMAT treatmentinvolving the first energy mode and the second energy mode.

In some embodiments, a trajectory optimizer may have a state that is afunction of the already discovered optimal trajectories. For example,for a given breast patient, a first optimal trajectory may have a firstcontrol point sequence S1 and energy 6×, a second optimal trajectory mayhave a second control point sequence S2 and energy 10× that complementsthe first optimal trajectory, and a third optimal trajectory with athird control point sequence S3 and energy 6× that complements the firstoptimal trajectory and the second optimal trajectory.

In some embodiments, considering M energy modes, the dose, score, andregion sections in the trajectory optimizer may be constructed such thateach connected component of the search graph corresponds to a singleenergy. As such, each trajectory will correspond to a fixed energy, andenergy switches will take place when transitioning from one trajectoryto another trajectory.

D. Second Method of Trajectory Optimization Considering MultipleCandidate Energy Modes

FIG. 15 is a flowchart illustrating a method 1500 of trajectoryoptimization considering multiple candidate energy modes according tosome other embodiments. The method 1500 may be referred to as amixed-energy-mode trajectory optimization over multiple candidate energymodes.

At 1502, a patient model is provided. The patient model includes one ormore regions of interest (ROIs) for the radiotherapy treatment.

At 1504, a delivery coordinate space (DCS) is defined. The DCS has a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane.

At 1506, a first energy mode and a second energy mode are identified forthe radiotherapy treatment.

At 1508, for the first energy mode, a first BEV region connectivitymanifold is determined by evaluating doses of the one or more ROIs foreach respective BEV plane corresponding to each respective candidatevertex of the DCS using transport solutions for the first energy mode.

At 1510, for the second energy mode, a second BEV region connectivitymanifold is determined by evaluating doses of the one or more ROIs foreach respective BEV plane corresponding to each respective candidatevertex of the DCS using transport solutions for the second energy mode.

At 1512, a set of optimal treatment trajectories is determined based onthe first BEV region connectivity manifold and the second BEV regionconnectivity manifold.

It should be appreciated that the specific steps illustrated in FIG. 15provide a particular method of trajectory optimization according to someembodiments. Other sequences of steps may also be performed according toalternative embodiments. For example, alternative embodiments mayperform the steps in a different order. Moreover, the individual stepsmay include multiple sub-steps that may be performed in varioussequences as appropriate to the individual step. Furthermore, additionalsteps may be added and some steps may be removed depending on theparticular applications. One of ordinary skill in the art wouldrecognize many variations, modifications, and alternatives.

E. Trajectory Optimization Including Consideration of Beam-Off Time

In the case of mixed-energy mode VMAT with N trajectories, the aggregatebeam-off time Δt=Σ_(i=1) ^(N−1) Δt_(i,i+1) may include, but not limitedto the following:

-   -   a time of Δt_(i,i+1) ^(switch) to switch between energy modes        m_(i) and m₁₊₁ and between trajectories i and i+1;    -   a time of Δt_(i,i+1) ^(move) to move machine axes from positions        at the last vertex v_(i) ^(last) of trajectory i to those at the        first vertex v_(i+1) ^(first) of trajectory i+1;    -   a time of Δt_(i,i+1) ^(collimator) to rotate the collimator from        its position c_(i) ^(last) at the last vertex of trajectory i to        its position c_(i+1) ^(first) at the first vertex of trajectory        i+1;    -   a time of Δt_(i,i+1) ^(imaging) of acquiring an image or        multiple images for guidance of treatment between the end of        trajectory i and the beginning of trajectory i+1.

The first item on the list arises from the introduction of multipleenergy modes.

There may exist a multitude of ways to minimize the aggregate beam-offtime Δt. In some embodiments, constraints may be imposed in amulti-energy mode VMAT optimization algorithm to minimize the temporaldistance between the last node of trajectory i, n_(last) ^((i)) and thefirst node of trajectory i+1, n_(first) ^((i+1)) minimal. This may beachieved by minimizing one or several of the constituents of Δt_(i,i+1)for each i∈{1, 2, . . . , N}. This may be referred to as theoptimization of the radiation order. In some other embodiments, theradiation order of the beams may be determined as a post-processing stepwithout the imposition of constraints between first and last nodes. Insome further embodiments, a combination of minimization internal to thealgorithm and post-processing may be used.

VIII. Beam Angle Optimization Considering Multiple Candidate EnergyModes

In some embodiments, the beam angle optimization (BAO) method describedabove may be generalized to consider multiple candidate energy modesusing two different approaches. In a first approach, a beam angleoptimization is performed for each of a plurality of candidate energymodes to obtain a respective candidate set of field geometries for eachcandidate energy mode. The candidate set of field geometries that givesrise to the optimal value of an objective function (“score”) may bepicked as the optimal set of field geometries. In some embodiments, theobjective function is the max-distance function in beam angleoptimization, as discussed above. The optimal set of field geometrieswill be used in the IMRT radiotherapy treatment at the energy modecorresponding to the optimal set of field geometries. This approach maybe referred to as single-energy-mode beam angle optimization overmultiple candidate energy modes.

In a second approach, the end result of a beam angle optimizationconsidering multiple candidate energy modes may be a set of fieldgeometries, possibly with different energy modes (e.g., some fieldgeometries in the set may correspond to a first energy mode, while someother field geometries in the set may correspond to a second energy modedifferent from the first energy mode). The second approach may bereferred to as mixed-energy-mode beam angle optimization over multiplecandidate energy modes.

A. Single-Energy-Mode Beam Angle Optimization Over Multiple CandidateEnergy Modes

According to some embodiments, to find an optimal energy mode for a setof k beams among M candidate energy modes (indexed as m=1, . . . M), theBAO method described above may be run M times, one for each of the Mcandidate energy modes. This process may produce M candidate sets of kbeams, each candidate set corresponding to a respective energy mode. Thecandidate set that gives rise to the optimal value of an objectivefunction (“score”) may be picked as the optimal set of k beams, whichwill be used for the IMRT radiotherapy treatment at the energy modecorresponding to the optimal set of k beams.

B. First Method of Beam Angle Optimization Considering MultipleCandidate Energy Modes

FIG. 16 is a flowchart illustrating a method 1600 of beam angleoptimization for an IMRT radiotherapy treatment considering multiplecandidate energy modes according to some other embodiments. The method1600 may be referred to as a sing-energy-mode beam angle optimizationover multiple candidate energy modes, as discussed above.

At 1602, a patient model is provided. The patient model includes one ormore regions of interest (ROIs) for the IMRT radiotherapy treatment.

At 1604, a delivery coordinate space (DCS) is defined. The DCS has a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane.

At 1606, a plurality of candidate energy modes is identified for theIMRT radiotherapy treatment.

At 1608, for each respective energy mode of the plurality of candidateenergy modes, for each respective BEV plane of a respective candidatevertex in the DCS, and for each respective ROI of the one or more ROIs:a respective dose of the respective ROI from a respective beamletincident from the respective candidate vertex and passing through arespective pixel of the respective BEV plane is evaluated usingtransport solutions for the respective energy mode; a respective BEVscore of the respective pixel is evaluated using the doses of the one ormore ROIs evaluated for the respective beamlet; one or more BEV regionsin the respective BEV plane for the respective energy mode aredetermined based on the BEV scores of the pixels of the BEV planescorresponding to the set of candidate vertices in the DCS; and arespective BEV region connectivity manifold for the respective energymode is determined based on the BEV regions of the BEV planes of the setof candidate vertices in the DCS. The respective BEV region connectivitymanifold represents connections between contiguous BEV regions betweenadjacent vertices.

At 1610, a plurality of candidate sets of IMRT fields is determined by:for each respective energy mode of the plurality of candidate energymodes, determining a respective candidate set of IMRT fields based onthe respective BEV region connectivity manifold for the respectiveenergy mode. Each respective IMRT field defines a beam anglecorresponding to a respective vertex in the DCS.

At 1612, one of the plurality of candidate sets of IMRT fields areselected as an optimal set of IMRT fields based on an objectivefunction. The optimal set of IMRT fields corresponds to an optimalenergy mode among the plurality of candidate energy modes.

It should be appreciated that the specific steps illustrated in FIG. 16provide a particular method of beam angle optimization for an IMRTradiotherapy treatment according to some embodiments. Other sequences ofsteps may also be performed according to alternative embodiments. Forexample, alternative embodiments may perform the steps in a differentorder. Moreover, the individual steps may include multiple sub-stepsthat may be performed in various sequences as appropriate to theindividual step. Furthermore, additional steps may be added and somesteps may be removed depending on the particular applications. One ofordinary skill in the art would recognize many variations,modifications, and alternatives.

C. Mixed-Energy-Mode Beam Angle Optimization Over Multiple CandidateEnergy Modes

According to some embodiments, to find a set of k beams among Mcandidate energy modes (indexed as m=1, . . . M), in which the k beamsmay have different energy modes, a multi-energy-mode BAO graph based onthe TORUS graph concept may be built based on the following nodes andedges:

Nodes

Nodes are a set

of N tuples of the form (energy mode, vertex, collimator index, regionmask). Each tuplet may determine a single field. An energy mode maycorrespond to a combination of photon energy and primary fluence mode. Avertex may correspond to a location in a discretized delivery coordinatespace (DCS), equivalent to specifying the isocenter, and gantry andcouch angles for a C-arm linear accelerator. A collimator index maycorrespond to a collimator angle out of a set of discrete possiblecollimator angles. A region mask may correspond a set of contiguoustarget regions the MLC leaves may expose.

Edges

Edges are connections between neighbor nodes that have vertex-vertexconnectivity in the underlying delivery coordinate space, and MLCconnectivity between the respective MLC leaf sequences. In someembodiments, there may be a connection between node n₁=(m₁, v, c, b) andnode n₂=(m₂, v, c, b) for any m₁ and m₂ among the M energy modes underconsideration. In some embodiments, there may be no edge between node(m₁, v₁, c₁, b₁) and node (m₂, v₂, c₂, b₂) if v₁≠v₂ or c₁≠c₂ or b₁≠b₂.In the case of beam-angle optimization, the MLC connectivity constraintmay not be important, but it may nevertheless act to reduce the numberof edges in the graph. Thus, it may be computationally useful to keepthe MLC connectivity constraint.

In the TORUS framework described above, min-distances and max-distancesare defined along edges in the TORUS graph. In the case of staticfields, there is no gantry motion while the treatment beam is on.Therefore, only the vertices themselves matter, and there is no need fora min-distance function. The score (max-distance) may be the metric tooptimize. The score S may be defined on subsets of nodes as,

S(

)=max−distance corresponding to the set of nodes (beams)

⊂

.  (18)

The optimal set of k beams may be defined to be the subset

⊂

with |

|=k that gives an optimal score S(

), where |

| represents the number of beams in the subset

⊂

.

The score may be a complex non-local function of the entire subset

. Therefore, it may be a non-trivial problem to solve exactly. However,there may exist an efficient way to find approximate solutions fasterthan trying all _(N)C_(k) combinations.

In some embodiments, a mixed-energy-mode-IMRT-plan BAO optimization mayinclude the following steps:

-   -   1. Evaluate the scores for individual beams (|        |=1), and use the squares of these scores as sampling        probabilities.    -   2. Randomly sample k beams using above probabilities.    -   3. Apply local gradient descent to this set        to find a local minimum.    -   4. Repeat steps 1-3 until there is no improvement in score found        for a predetermined number of successive trials (e.g., 20        successive trials).

This procedure may be referred to as a gradient descent method. Duringeach iteration of the gradient descent procedure, all possible “edges”may be considered looking for improvement, where an “edge” is motion ofthe subset

to a neighbor subset

′. The subset of beams may be updated to move in the direction that mostimproves the score. This is repeated until no more local improvements inthe score is found. A neighbor to subset

is defined to be another subset B′ that differs in exactly one beamB_(i)→B′_(i), where these two beams either share an edge in the TORUSgraph or have the same position in the delivery coordinate space (samevertex).

D. Second Method of Beam Angle Optimization Considering MultipleCandidate Energy Modes

FIG. 17 is a flowchart illustrating a method 1700 of beam angleoptimization for an IMRT radiotherapy treatment considering multiplecandidate energy modes according to some other embodiments. The method1700 may be referred to as mixed-energy-mode beam angle optimizationover multiple candidate energy modes, as discussed above.

At 1702, a patient model is provided. The patient model includes one ormore regions of interest (ROIs) for the IMRT radiotherapy treatment.

At 1704, a delivery coordinate space (DCS) is defined. The DCS has a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane.

At 1706, a first energy mode and a second energy mode are identified forthe IMRT radiotherapy treatment.

At 1708, for the first energy mode, a first BEV region connectivitymanifold is determined by evaluating doses of the one or more ROIs foreach respective BEV plane corresponding to each respective candidatevertex of the DCS using transport solutions for the first energy mode.

At 1710, for the second energy mode, a second BEV region connectivitymanifold is determined by evaluating doses of the one or more ROIs foreach respective BEV plane corresponding to each respective candidatevertex of the DCS using transport solutions for the second energy mode.

At 1720, a set of IMRT fields is determined based on the first BEVregion connectivity manifold for the first energy mode and the secondBEV region connectivity manifold for the second energy mode. Eachrespective IMRT field of the set of IMRT fields defines a beam anglecorresponding to a respective vertex in the DCS.

It should be appreciated that the specific steps illustrated in FIG. 17provide a particular method of beam angle optimization for an IMRTradiotherapy treatment according to some embodiments. Other sequences ofsteps may also be performed according to alternative embodiments. Forexample, alternative embodiments may perform the steps in a differentorder. Moreover, the individual steps may include multiple sub-stepsthat may be performed in various sequences as appropriate to theindividual step. Furthermore, additional steps may be added and somesteps may be removed depending on the particular applications. One ofordinary skill in the art would recognize many variations,modifications, and alternatives.

E. Beam angle optimization including consideration of beam-off time

Consider a multi-energy mode IMRT treatment as an external-beamradiotherapy treatment, where the dose is delivered from k static beamlocations r_(i)=(m_(i), v_(i), c_(i)), i=1, . . . , k, where the indicesm_(i), v_(i), and c_(i) correspond to the energy mode, vertex, andcollimator angle, respectively. In a C-arm linear accelerator treatmentsystem, a vertex may include an isocenter, a gantry angle, and a couchangle. In other types of external-beam radiation treatment systems, avertex may include other treatment axes variables. The fluence deliveredfrom each r_(i) may be determined by way of a fluence-optimizationscheme that distributes the fluence optimally between the r_(i) suchthat clinical optimization objectives are fulfilled.

A dosimetrically optimal IMRT treatment plan with multi-energy modes maybe found within a given delivery coordinate space by running fluenceoptimization and leaf sequencing for each permissible combination of{r_(i)}_(i=1) ^(k), and by picking the dosimetrically optimal one.However, this process may be impractical due to the large number ofpossible combinations. Moreover, it may be desirable that the patientand the patient's internal organs remain stationary during theadministration of radiation for the delivered dose to match the planneddose. The longer the treatment takes, the more likely that the patientor the patient's internal organs may move during the treatment, andhence a higher probability of not delivering the intended dose to thetarget volumes.

The overall treatment time may be prolonged by the beam-off transitiontimes from beam r_(i) to r_(i+1). The total beam-off transition timefrom beam r_(i) to r_(i+1) may be expressed asΔt_(i,i+1)=max_(j)Δt_(i,i+1) ^(change) ^(j) , where the changes mayinclude, but are limited to the following:

-   -   a time of Δt_(i,i+1) ^(switch) to switch between energy modes        m_(i) and m_(i+1);    -   a time of Δt_(i,i+1) ^(move) to move machine axes from positions        at vertex v_(i) to those at v_(i+1);    -   a time of Δt_(i,i+1) ^(collimator) to rotate the collimator from        c_(i) to c_(i+1);    -   a time of Δt_(i,i+1) ^(imaging) of acquiring an image or        multiple images for guidance of treatment between beams i and        i+1.

The first item on the list arises from the introduction of multipleenergy modes.

A minimization of the aggregate beam-off time Δt=Σ_(i=1) ^(k−1)Δt_(i,i+1) may be an important part of the beam-selection problem. Insome embodiments, a multi-energy BAO algorithm may include constraintson picking beams that can be delivered in a time-efficient radiationorder. In some other embodiments, the radiation order of the beams maybe determined as a post-processing step.

IX. Hybrid Trajectory and Beam Angle Optimization

The TORUS approach described above may find continuous trajectoriesamenable to VMAT-like delivery in a given delivery coordinate space. Therequirement of spatial continuity in TORUS may come at the potentialexpense of reduction in dosimetric plan quality. To overcome thispotential drawback, the VMAT trajectories may be augmented by placingone or more IMRT fields, either in the same delivery coordinate space asthe VMAT trajectories or in a different delivery coordinate space. TheIMRT fields may be obtained using the methods described above in SectionVI. The resulting treatment geometry may be referred to as a hybridVMAT-IMRT treatment geometry.

Hybrid VMAT-IMRT treatment geometries may be generated in a multitude ofways. In the following, cases that are either inherently time-efficientor may provide a significant improvement in dosimetric plan quality atthe expense of longer treatment times according to some embodiments arediscussed.

A. VMAT Treatment Geometry Augmented with IMRT Fields in the SameDelivery Coordinate Space

Assume that a VMAT treatment geometry of N trajectories has beengenerated in a given delivery coordinate space, which may be abbreviatedas DCS₁, using any of the procedures described above in Section V for agiven energy mode. A given delivery coordinate space includes a set ofvertices that represent permissible directions of incidence that wouldnot result in collisions (e.g., machine-to-machine collisions andmachine-to-patient collisions). In some embodiments, the VMAT treatmentgeometry may be augmented by placing k IMRT fields along one or moretrajectories of the N trajectories. The k IMRT fields may be obtainedusing any of the methods described above in Section VI. In some otherembodiments, the angles of the k IMRT fields along the VMAT-liketrajectories may be selected by means different from those describedabove in Section VI, in conjunction with concurrent optimization of MUtime series and leaf positioning.

As an example, consider a C-arm linear accelerator radiotherapytreatment system illustrated in FIG. 18A. The radiotherapy treatmentsystem includes a couch 1802 for supporting a patient 1806 (a patientdummy is shown), and a gantry 1804 that rotates around the couch 1802.Consider a first delivery coordinate space DCS₁, which has a fixedisocenter, a couch angle of 0 deg, and a gantry angle interval of [0deg, 360 deg), as illustrated in FIG. 18A. DCS₁ may be referred to as a“coplanar” delivery coordinate space. FIG. 18B shows a view of a patientmodel 1810 from the top of the patient's head and a VMAT-like arc 1820.In this case, the VMAT-like arc 1820 happens to a complete circle aroundthe patient's head (e.g., for treating intracranial metastasis). TheVMAT-like arc 1820 may be obtained by a TORUS trajectory optimizationdescribed above in Section V.

In some embodiments, one or more IMRT fields may be inserted at certaindirections along the VMAT-like arc 1820 to augment the VMAT-like arc1820. The directions of the IMRT fields may be determined using a beamangle optimization described above in Section VI. Alternatively, aradiologist may determine that delivering more radiation doses at thosedirections may improve the dosimetry quality of the radiotherapytreatment. Such a hybrid VMAT-IMRT treatment geometry may be efficientin terms of total delivery time. For example, the radiotherapy treatmentsystem may start delivery of the VMAT-like arc 1820 by rotating thegantry 18004 along the VMAT-like arc 1820, then stop at the IMRT angleto deliver the IMRT field, and then resume rotating the gantry 1804 fromthat angle to complete the VMAT-like arc 1820. In this manner, theradiotherapy treatment may not incur any beam-off transition timebetween delivery of the VMAT-like arc 1820 and the delivery of the IMRTfield.

In some embodiments, a TORUS trajectory optimization may generateadditional VMAT-like arcs in the same delivery coordinate space DCS₁.For example, as shown in FIG. 18B, a second VMAT-like arc 1830 may begenerated. In this case, the second VMAT-like arc 1830 is not a completecircle, but is missing a section in the fourth quadrant on the upperleft. (The second VMAT-like arc 1830 may have an MLC leaf sequence whoseapertures may overlap fully, partially or not at all with those of thefirst VMAT-like arc 1820 at vertices common to both VMAT-like arcs 1820and 1830.)

It is possible that a beam angle optimization may identify one or moreoptimal directions of IMRT fields (e.g., the direction indicated by thedashed line 1840 in FIG. 18B) that fall in the missing quadrant (notethose angles are still in the same delivery coordinate space DCS₁).

In some embodiments, if a VMAT treatment geometry includes a first arc1820 and a second arc 1830, and the direction 1840 of the IMRT fieldfalls on only the first arc 1820, the IMRT field may be delivered inconjunction with the delivery of the first arc 1820, instead of inconjunction with the delivery of the second arc 1830. This ordering ofthe hybrid VMAT-IMRT treatment geometry may be more time-efficient thanthe case if the IMRT field is delivered in conjunction with the deliveryof the second arc 1830. This is because, in the latter case, theradiotherapy treatment system may need to turn off the radiation beam(either in the middle of the second arc 1830 or after completing thesecond arc 1830), rotate the gantry 1804 to the direction 1840 of theIMRT field, then deliver the IMRT field. Therefore, the radiotherapytreatment may incur some beam-off transition time.

In some embodiments, it may happen that an optimal direction of a IMRTfield does not fall on any of the N trajectories of the VMAT treatmentgeometry. For instance, in the example illustrated in FIG. 18C, assumethat the VMAT treatment geometry includes two VMAT-like arcs 1850 and1860. The combination of the two arcs 1850 and 1860 does not cover acomplete circle. It may be possible that an optimal direction 1870 ofthe IMRT field may fall on neither of the two arcs 1850 and 1860. Insome embodiments, such a hybrid VMAT-IMRT treatment geometry may stillbe preferred, if the IMRT field provides a significant improvement indosimetric plan quality, or it may improve the patient's quality of life(e.g., if the patient's eyesight may be saved), albeit at the expense oflonger treatment time.

In some embodiments, multiple treatment axes may rotate and/or translateconcurrently in a VMAT-like trajectory. The treatment axes may include,for example, the couch position (in three degrees of freedom), the couchangle (pitch, yaw, and roll), the collimator angle, and the like, inaddition to the gantry angle. For instance, the couch can be rotatedsimultaneously with the rotation of the gantry. The collimator may bealso be rotated at the same time. In some embodiments, a TORUStrajectory optimization algorithm may generate such multi-axestrajectories.

FIG. 19 shows an example of a more general VMAT-like trajectory 1920superimposed on a patient model 1910. Similar to the examples discussedabove in conjunction with FIGS. 18B, one or more IMRT fields may beinserted along certain directions along the VMAT-like trajectory 1920,to result in a time-efficient hybrid VMAT-IMRT treatment geometry.Alternatively, a hybrid VMAT-IMRT treatment geometry may include IMRTfields that do not fall on the VMAT-like trajectory 1920, if those IMRTfields may improve the dosimetric plan quality or improve the patient'squality of life.

It should be understood that the methods discussed above are not limitedto the couch model illustrated in FIG. 18A, but can be applied to othercouch models. The methods are also not limited to C-arm linearaccelerator radiotherapy treatment systems, but can be applied to othertypes of radiotherapy treatment systems.

B. First Method of Hybrid VMAT-IMRT Treatment Geometry Optimization

FIG. 20 is a flowchart illustrating a method 2000 of hybrid VMAT-IMRTtreatment geometry optimization according to some embodiments.

At 2002, a patient model is provided. The patient model includes one ormore regions of interest (ROIs) for the radiotherapy treatment.

At 2004, a delivery coordinate space (DCS) is defined. The DCS has a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane.

At 2006, for each respective BEV plane of a respective candidate vertexin the DCS, and for each respective ROI of the one or more ROIs: arespective dose of the respective ROI from a respective beamlet incidentfrom the respective candidate vertex and passing through a respectivepixel of the respective BEV plane is evaluated using transport solutionsof the respective beamlet; a respective BEV score of the respectivepixel is evaluated using the doses of the one or more ROIs evaluated forthe respective beamlet; and one or more BEV regions in the respectiveBEV plane are determining based on the BEV scores of the pixels of theBEV planes corresponding to the set of candidate vertices in the DCS.

At 2008, a BEV region connectivity manifold is determined based on theBEV regions of the BEV planes of the set of candidate vertices in theDCS. The BEV region connectivity manifold represents connections betweencontiguous BEV regions between adjacent vertices.

At 2010, a set of treatment trajectories is determined based on the BEVregion connectivity manifold. Each treatment trajectory defines arespective path through a respective set of vertices in the DCS.

At 2012, one or more IMRT fields are determined. Each respective IMRTfield defines a respective direction of incidence corresponding to arespective vertex in the DCS.

It should be appreciated that the specific steps illustrated in FIG. 20provide a particular method of hybrid VMAT-IMRT treatment geometryoptimization according to some embodiments. Other sequences of steps mayalso be performed according to alternative embodiments. For example,alternative embodiments may perform the steps in a different order.Moreover, the individual steps may include multiple sub-steps that maybe performed in various sequences as appropriate to the individual step.Furthermore, additional steps may be added and some steps may be removeddepending on the particular applications. One of ordinary skill in theart would recognize many variations, modifications, and alternatives.

C. VMAT Treatment Geometry Augmented with IMRT Fields on a DifferentDelivery Coordinate Space

In some embodiments, VMAT-like trajectories may be augmented by placingk IMRT fields generated using the methods described in Section VI in theunion of DCS₁ and another delivery coordinate space, DCS₂. As anexample, consider a C-arm linear accelerator. A first deliverycoordinate space DCS₁ may have a fixed isocenter, a couch angle of 0deg, and a gantry angle interval of [0 deg, 360 deg), as illustrated inFIG. 21A. A second delivery coordinate space DCS₂ may have the sameisocenter as DCS₁, a couch angle of 90 deg with an associated gantryangle interval of (180 deg, 360 deg], or a couch angle of 270 deg withan associated gantry angle interval of [0 deg, 180 deg), as illustratedin FIG. 21C. DCS₁ and DCS₂ are two coplanar delivery coordinate spacesbut on different planes, and the only vertex they share has a gantryangle of 0 deg since rotating the couch does not change the direction ofincidence when the gantry angle is equal to 0 deg. FIG. 21B shows aVMAT-like trajectory 2120 in DCS₁ (which happens to be a completecircle) superimposed on a patient model 2110. DCS₂ may provide anon-coplanar extension for shaping the entry and exit doses about thetreatment target. Thus, in some embodiments, one or more IMRT fields maybe added in DCS₂ to augment the VMAT treatment geometry in DCS₁, forminga hybrid VMAT-IMRT treatment geometry in two different deliverycoordinate spaces.

In some embodiments, a VMAT-IMRT treatment geometry may includeadditional VMAT-like trajectories in additional delivery coordinatespaces. For example, FIG. 21D shows a VMAT-like trajectory 2130 in DCS₂.FIG. 22 shows two VMAT-like trajectories 2220 and 2230 in a thirddelivery coordinate space DCS₃ and a fourth delivery coordinate spaceDCS₄, respectively. Thus, in some exemplary embodiments, a VMATtreatment geometry may include four (or more) VMAT-like trajectories infour (or more) different delivery coordinate spaces. While it may bemore complicated to deliver a radiotherapy treatment in multipledelivery coordinate spaces, such multi-DCS treatment geometry may beselected if it may improve the dosimetric plan quality.

D. Second Method of Hybrid VMAT-IMRT Treatment Geometry Optimization

FIG. 23 is a flowchart illustrating a method 2300 of hybrid VMAT-IMRTtreatment geometry optimization according to some embodiments.

At 2302, a patient model is provided. The patient model includes one ormore regions of interest (ROIs) for the radiotherapy treatment.

At 2304, a first delivery coordinate space and a second deliverycoordinate space is defined. The first delivery coordinate space has afirst set of candidate vertices. The second delivery coordinate spacehas a second set of candidate vertices. Each vertex of the first set ofcandidate vertices or the second set of candidate vertices defines arespective beam's eye view (BEV) plane.

At 2306, a first beam's eye view (BEV) region connectivity manifold isdetermined by evaluating doses of the one or more ROIs for eachrespective BEV plane corresponding to each respective vertex of thefirst set of candidate vertices of the first delivery coordinate space.

At 2308, a first set of treatment trajectories is determined based onthe first BEV region connectivity manifold. Each treatment trajectory ofthe first set of treatment trajectories defines a respective paththrough a respective set of vertices in the first delivery coordinatespace.

At 2310, a first set of IMRT fields is determined. Each of the first setof IMRT fields corresponds to a respective vertex in the second deliverycoordinate space.

It should be appreciated that the specific steps illustrated in FIG. 23provide a particular method of hybrid VMAT-IMRT treatment geometryoptimization according to some embodiments. Other sequences of steps mayalso be performed according to alternative embodiments. For example,alternative embodiments may perform the steps in a different order.Moreover, the individual steps may include multiple sub-steps that maybe performed in various sequences as appropriate to the individual step.Furthermore, additional steps may be added and some steps may be removeddepending on the particular applications. One of ordinary skill in theart would recognize many variations, modifications, and alternatives.

E. VMAT Treatment Geometry Augmented with IMRT Fields on the SameDelivery Coordinate Space with Mixed Energy Modes

Assume that a VMAT-treatment geometry of N trajectories includetrajectories of different energy modes. For example, as discussed abovein Section VII.C., a trajectory optimization algorithm may generate sometrajectories of a first energy mode and some other trajectories in asecond energy mode different from the first energy mode. In someembodiments, the VMAT-treatment geometry may be augmented by placing kIMRT fields along the trajectories, such that the energy mode of eachIMRT field corresponds to the energy mode of the VMAT-like trajectory onwhich the IMRT field sits. The k IMRT fields may be obtained using anyof the methods described above in Sections VI and VIII. This method maynot incur additional treatment delivery time due to energy switchingbetween VMAT and IMRT fields. Note also that the positions of the IMRTfields may be set prior to concurrent optimization of MU time series andleaf positioning.

In some other embodiments, the positions of the k IMRT fields along theVMAT-like trajectories may be selected by means different from thosedescribed above in Section VI, in conjunction with concurrentoptimization of MU time series and leaf positioning.

In some embodiments, within DCS₁, but outside of trajectories, the setof viable energy modes for the IMRT fields may be determined based onother constraints such as total delivery time of treatment.

F. Third Method of Hybrid VMAT-IMRT Treatment Geometry Optimization

FIG. 24 is a flowchart illustrating a method 2400 of hybrid VMAT-IMRTtreatment geometry optimization according to some embodiments.

At 2402, a patient model is provided. The patient model includes one ormore regions of interest (ROIs) for the radiotherapy treatment.

At 2404, a delivery coordinate space (DCS) is defined. The DCS has a setof candidate vertices. Each respective candidate vertex defines arespective beam's eye view (BEV) plane.

At 2406, a first energy mode and a second energy mode are identified forthe radiotherapy treatment.

At 2408, for the first energy mode, determine a first BEV regionconnectivity manifold by evaluating doses of the one or more ROIs foreach respective BEV plane corresponding to each respective candidatevertex of the DCS using transport solutions for the first energy mode.

At 2410, for the second energy mode, determine a second BEV regionconnectivity manifold by evaluating doses of the one or more ROIs foreach respective BEV plane corresponding to each respective candidatevertex of the DCS using transport solutions for the second energy mode.

At 2412, a set of optimal treatment trajectories is determined based onthe first BEV region connectivity manifold and the second BEV regionconnectivity manifold.

At 2414, a set of IMRT fields is determined. Each respective IMRT fieldof the set of IMRT fields corresponds to a respective vertex in thedelivery coordinate space.

It should be appreciated that the specific steps illustrated in FIG. 24provide a particular method of hybrid VMAT-IMRT treatment geometryoptimization according to some embodiments. Other sequences of steps mayalso be performed according to alternative embodiments. For example,alternative embodiments may perform the steps in a different order.Moreover, the individual steps may include multiple sub-steps that maybe performed in various sequences as appropriate to the individual step.Furthermore, additional steps may be added and some steps may be removeddepending on the particular applications. One of ordinary skill in theart would recognize many variations, modifications, and alternatives.

X. Computer System

Any of the computer systems mentioned herein may utilize any suitablenumber of subsystems. Examples of such subsystems are shown in FIG. 25in computer system 2000. In some embodiments, a computer system includesa single computer apparatus, where the subsystems can be the componentsof the computer apparatus. In other embodiments, a computer system caninclude multiple computer apparatuses, each being a subsystem, withinternal components.

The subsystems shown in FIG. 25 are interconnected via a system bus2575. Additional subsystems such as a printer 2574, keyboard 2578,storage device(s) 2579, monitor 2576, which is coupled to displayadapter 2582, and others are shown. Peripherals and input/output (I/O)devices, which couple to I/O controller 2571, can be connected to thecomputer system by any number of means known in the art, such as serialport 2577. For example, serial port 2577 or external interface 2581(e.g. Ethernet, Wi-Fi, etc.) can be used to connect computer system 2500to a wide area network such as the Internet, a mouse input device, or ascanner. The interconnection via system bus 2575 allows the centralprocessor 2573 to communicate with each subsystem and to control theexecution of instructions from system memory 2572 or the storagedevice(s) 2579 (e.g., a fixed disk, such as a hard drive or opticaldisk), as well as the exchange of information between subsystems. Thesystem memory 2572 and/or the storage device(s) 2579 may embody acomputer readable medium. Any of the data mentioned herein can be outputfrom one component to another component and can be output to the user.

A computer system can include a plurality of the same components orsubsystems, e.g., connected together by external interface 2581 or by aninternal interface. In some embodiments, computer systems, subsystem, orapparatuses can communicate over a network. In such instances, onecomputer can be considered a client and another computer a server, whereeach can be part of a same computer system. A client and a server caneach include multiple systems, subsystems, or components.

External interface 2581 can be used to transmit one or more treatmentplans to one or more radiation treatment devices, as described herein.For example, a treatment planning application can reside on a servercomputer, and a client computer can use the treatment planningapplication. The server computer can be part of a cloud computingplatform that provides software as a service (SaaS). Once a treatmentplan is determined, a client computer can specify which radiation deviceor a treatment plan database accessible by the radiation device fortransmitting one or more files encapsulating the treatment plan. Forinstance, an IP address can be specified.

It should be understood that any of the embodiments of the presentinvention can be implemented in the form of control logic using hardware(e.g. an application specific integrated circuit or field programmablegate array) and/or using computer software with a generally programmableprocessor in a modular or integrated manner. As used herein, a processorincludes a multi-core processor on a same integrated chip, or multipleprocessing units on a single circuit board or networked. Based on thedisclosure and teachings provided herein, a person of ordinary skill inthe art will know and appreciate other ways and/or methods to implementembodiments of the present invention using hardware and a combination ofhardware and software.

Any of the software components or functions described in thisapplication may be implemented as software code to be executed by aprocessor using any suitable computer language such as, for example,Java, C++ or Perl using, for example, conventional or object-orientedtechniques. The software code may be stored as a series of instructionsor commands on a computer readable medium for storage and/ortransmission, suitable media include random access memory (RAM), a readonly memory (ROM), a magnetic medium such as a hard-drive or a floppydisk, or an optical medium such as a compact disk (CD) or DVD (digitalversatile disk), flash memory, and the like. The computer readablemedium may be any combination of such storage or transmission devices.

Such programs may also be encoded and transmitted using carrier signalsadapted for transmission via wired, optical, and/or wireless networksconforming to a variety of protocols, including the Internet. As such, acomputer readable medium according to an embodiment of the presentinvention may be created using a data signal encoded with such programs.Computer readable media encoded with the program code may be packagedwith a compatible device or provided separately from other devices(e.g., via Internet download). Any such computer readable medium mayreside on or within a single computer product (e.g. a hard drive, a CD,or an entire computer system), and may be present on or within differentcomputer products within a system or network. A computer system mayinclude a monitor, printer, or other suitable display for providing anyof the results mentioned herein to a user.

Any of the methods described herein may be totally or partiallyperformed with a computer system including one or more processors, whichcan be configured to perform the steps. Thus, embodiments can bedirected to computer systems configured to perform the steps of any ofthe methods described herein, potentially with different componentsperforming a respective steps or a respective group of steps. Althoughpresented as numbered steps, steps of methods herein can be performed ata same time or in a different order. Additionally, portions of thesesteps may be used with portions of other steps from other methods. Also,all or portions of a step may be optional. Additionally, any of thesteps of any of the methods can be performed with modules, circuits, orother means for performing these steps.

The specific details of particular embodiments may be combined in anysuitable manner without departing from the spirit and scope ofembodiments of the invention. However, other embodiments of theinvention may be directed to specific embodiments relating to eachindividual aspect, or specific combinations of these individual aspects.

The above description of exemplary embodiments of the invention has beenpresented for the purposes of illustration and description. It is notintended to be exhaustive or to limit the invention to the precise formdescribed, and many modifications and variations are possible in lightof the teaching above. The embodiments were chosen and described inorder to best explain the principles of the invention and its practicalapplications to thereby enable others skilled in the art to best utilizethe invention in various embodiments and with various modifications asare suited to the particular use contemplated.

A recitation of “a”, “an” or “the” is intended to mean “one or more”unless specifically indicated to the contrary.

All patents, patent applications, publications, and descriptionsmentioned herein are incorporated by reference in their entirety for allpurposes. None is admitted to be prior art.

1-9. (canceled)
 10. A method of beam angle optimization for an IMRTradiotherapy treatment, the method comprising: providing a patient modelincluding one or more regions of interest (ROIs) for the IMRTradiotherapy treatment; defining a delivery coordinate space (DCS)having a set of candidate vertices, each respective candidate vertexdefining a respective beam's eye view (BEV) plane; identifying aplurality of candidate energy modes for the IMRT radiotherapy treatment;for each respective energy mode of the plurality of candidate energymodes: for each respective BEV plane of a respective candidate vertex inthe DCS: for each respective ROI of the one or more ROIs: evaluating arespective dose of the respective ROI from a respective beamlet incidentfrom the respective candidate vertex and passing through a respectivepixel of the respective BEV plane using transport solutions for therespective energy mode; evaluating a respective BEV score of therespective pixel using the doses of the one or more ROIs evaluated forthe respective beamlet; and determining one or more BEV regions in therespective BEV plane for the respective energy mode based on the BEVscores of the pixels of the BEV planes corresponding to the set ofcandidate vertices in the DOS; and determining a respective BEV regionconnectivity manifold for the respective energy mode based on the BEVregions of the BEV planes of the set of candidate vertices in the DCS,the respective BEV region connectivity manifold representing connectionsbetween contiguous BEV regions between adjacent vertices; determining aplurality of candidate sets of IMRT fields by: for each respectiveenergy mode of the plurality of candidate energy modes: determining arespective candidate set of IMRT fields based on the respective BEVregion connectivity manifold for the respective energy mode, eachrespective IMRT field defining a beam angle corresponding to arespective vertex in the DOS; and selecting one of the plurality ofcandidate sets of IMRT fields as an optimal set of IMRT fields based onan objective function, the optimal set of IMRT fields corresponding toan optimal energy mode among the plurality of candidate energy modes.11. The method of claim 10, wherein the transport solutions for eachenergy mode is obtained by using an adjoint transport approach.
 12. Themethod of claim 10, wherein the transport solutions for each energy modeis obtained by using a forward transport approach.
 13. The method ofclaim 10, wherein determining the one or more BEV regions in therespective BEV plane for the respective energy mode comprises:determining a respective threshold BEV score for the respective energymode based on the BEV scores of the pixels of the BEV planescorresponding to the set of candidate vertices in the DOS; anddetermining the one or more BEV regions in the respective BEV plane bycomparing each respective BEV score of a respective pixel of therespective BEV plane to the respective threshold BEV score for therespective energy mode, wherein pixels within the one or more BEVregions have BEV scores greater than or equal to the respectivethreshold BEV score.
 14. The method of claim 10, wherein the pluralityof ROIs includes one or more planning target volumes (PTVs) and one ormore organs at risk (OARs), and evaluating the respective BEV score ofthe respective pixel comprises: evaluating a weighted linear combinationof the doses of the plurality of ROIs, wherein the dose of a respectiveROI that is one of the one or more PTVs has a positive weight, and adose of a respective ROI that is one of the one or more OARs has anegative weight.
 15. The method of claim 10, wherein determining therespective candidate set of IMRT fields comprises performingoptimization using a max-distance function using local gradient descentalgorithm based on information contained in the respective BEV regionconnectivity manifold for the respective energy mode.
 16. A method ofbeam angle optimization in an IMRT radiotherapy treatment, the methodcomprising: providing a patient model including one or more regions ofinterest (ROIs) for the IMRT radiotherapy treatment; defining a deliverycoordinate space (DCS) having a set of candidate vertices, eachrespective candidate vertex defining a respective beam's eye view (BEV)plane; identifying a first energy mode and a second energy mode for theIMRT radiotherapy treatment; for the first energy mode: determining afirst BEV region connectivity manifold by evaluating doses of the one ormore ROIs for each respective BEV plane corresponding to each respectivecandidate vertex of the DCS using transport solutions for the firstenergy mode; for the second energy mode: determining a second BEV regionconnectivity manifold by evaluating doses of the one or more ROIs foreach respective BEV plane corresponding to each respective candidatevertex of the DCS using transport solutions for the second energy mode;and determining a set of IMRT fields based on the first BEV regionconnectivity manifold for the first energy mode and the second BEVregion connectivity manifold for the second energy mode, each respectiveIMRT field of the set of IMRT fields defining a beam angle correspondingto a respective vertex in the DCS.
 17. The method of claim 16, whereinthe set of IMRT fields includes at least a first IMRT fieldcorresponding to the first energy mode, and at least a second IMRT fieldcorresponding to the second energy mode.
 18. The method of claim 16,wherein: determining the first BEV region connectivity manifold or thesecond BEV region connectivity manifold comprises: for each respectiveBEV plane of a respective candidate vertex in the DCS: for eachrespective ROI of the one or more ROIs: evaluating a respective dose ofthe respective ROI from a respective beamlet incident from therespective candidate vertex and passing through a respective pixel ofthe respective BEV plane using the transport solutions for the firstenergy mode or the second energy mode; and evaluating a respective BEVscore of the respective pixel using the doses of the one or more ROIsevaluated for the respective beamlet; and determining one or more BEVregions in the respective BEV plane for the first energy mode or thesecond energy mode based on the BEV scores of the pixels of the BEVplanes corresponding to the set of candidate vertices in the DOS; anddetermining the first BEV region connectivity manifold for the firstenergy mode based on the BEV regions for the first energy mode of theBEV planes of the set of candidate vertices in the DCS, or determiningthe second BEV region connectivity manifold for the second energy modebased on the BEV regions for the second energy mode of the BEV planes ofthe set of candidate vertices in the DCS.
 19. The method of claim 16,wherein the transport solutions for the first energy mode or the secondenergy mode is obtained by using an adjoint transport approach.
 20. Themethod of claim 16, wherein the transport solutions for the first energymode or the second energy mode is obtained by using a forward transportapproach.